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Abstract 

We present a simple two-dimensional dynamical system where two nonlinear terms, exerting respectively posi- 
tive feedback and reversal, compete to create a singularity in finite time decorated by accelerating oscillations. The 
power law singularity results from the increasing growth rate. The oscillations result from the restoring mechanism. 
As a function of the order of the nonlinearity of the growth rate and of the restoring term, a rich variety of behavior is 
documented analytically and numerically. The dynamical behavior is ttaced back fundamentally to the self-similar 
spiral structure of trajectories in phase space unfolding around an unstable spiral point at the origin. The interplay 
between the restoring mechanism and the nonlinear growth rate leads to approximately log-periodic oscillations 
with remarkable scaling properties. Three domains of applications are discussed: (1) the stock market with a com- 
petition between nonlinear ttend-followers and nonlinear value investors; (2) the world human population with a 
competition between a population-dependent growth rate and a nonlinear dependence on a finite carrying capacity; 
(3) the failure of a material subjected to a time-varying stress with a competition between positive geometrical 
feedback on the damage variable and nonlinear healing. 
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1 Introduction 



The mathematics of singularities is applied routinely in the physics of phase transitions to describe for instance the 
transformations from ice to water or from a magnet to a demagnetized state when raising the temperature, as well 
as in many other condensed matter systems. Such singularities characterize so-called critical phenomena. In these 
problems, physical observables such as susceptibilities, specific heat, etc., exhibit a singularity as the control parameter 
(temperature, strength of the interaction) approaches a critical value. 

Other classes of singularities occur in dynamical systems and are spontaneously reached in finite time. Sponta- 
neous singularities in ordinary (ODE) and partial differential equations (PDE) are quite common and have been found 
in many well-established models of natural systems, either at special points in space such as in the Euler equations 
of inviscid fluids p5| , |J, in the surface curvature on the free surface of a conducting fluid in an electric field [f59|], 
in vortex collapse of systems of point vortices, in the equations of General Relativity coupled to a mass field leading 
to the formation of black holes JTT|], in models of micro-organisms aggregating to form fruiting bodies [fl6|], or in 



the more prosaic rotating coin (Euler 's disk) [42]. Some more complex examples are models of rupture and material 
failure [25, 29], earthquakes [32] and stock market crashes [pUffiH]. 
The normal form of a finite-time singularity is the equation 



— = p m , with m > 1 , (1) 



dt 

whose solution is 

p(t)=p(0) (JL— J , (2) 

where the critical time t c = {m — l)/[p(0)] m_1 is determined by the initial condition p(0). The singularity results 
from the fact that the instantaneous growth rate dhxp/dt = p m ~ x is increasing with p and thus with time. This 
can be visualized by studying the doubling time, defined at the time interval At necessary for p(t) to double, i.e., 
p(t + At) = 2p(t). When the growth rate of p increases as a power law of p, the doubling time decreases fast and 
the sequence of doubling time intervals shrinks to zero sufficiently fast so that its sum is a convergent geometrical 
series. The variable thus undergoes an infinite number of doubling operations in a finite time, which the essence of a 
finite-time singularity. 

The power law solution (§) possesses the symmetry of "scale invariance", namely a reduction t c — t — > (t c — t) /A 
of the distance t c — t from the singularity at t c by an arbitrary factor A changes p(t) to \ l K' m ~ l ) p(t), i.e., keeps the 
same form of the solution up to a global rescaling. 

This continuous scale invariance can be partially broken into a weaker symmetry, called discrete scale invariance, 
according to which the self-similarity holds only for integer powers of a specific factor A p^Q. The hallmark of 
this discrete scale invariance is that the power law (||) transforms into an oscillatory singularity, with log-periodic 
oscillations decorating the overall power law acceleration. Such log -periodic power laws have been documented 
for many systems such as with a built-in geometrical hierarchy, in programming and number theory, for Newcomb- 
Benford law of first digits and in the arithmetic system, in diffusion in anisotropic quenched random lattices, as the 
result of a cascade of ultra-violet instabilities in growth processes and rupture, in deterministic dynamical systems 
(cascades of sub-harmonic bifurcations in the transition to chaos, two-coupled anharmonic oscillators, near-separatrix 
Hamiltonian chaotic dynamics, kicked charged particle moving in a double-well potential giving a physical realization 
of Mandelbrot and Julia sets, chaotic scattering), in extension of percolation theory (so-called "animals"), in response 
functions of spin systems with quenched disorder, in freely decaying 2D-turbulence, in the gravitational collapse and 



black hole formation, in spinodal decomposition of binary mixtures in uniform shear flow, etc. (see p9| , plTj and 
references therein). 

The novel interesting feature is the presence of a discrete hierarchy of length and/or time scales in an otherwise 
scale-invariant system. The presence of these scales may provide insight into the underlying mechanisms. While 
there is a good general framework for the description of discrete scale invariant systems using renormalization group 
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theory [|49|], a general understanding of the possible physical mechanisms at its origin is still lacking. In particular, 
dynamically generated discrete scale invariance is the most important problem, as it might provide understanding in 
the origins of the ubiquitous existence of hierarchies and cascades in natural and social systems. 

Here, we introduce and study a simple two-dimensional nonlinear dynamical system with the minimal ingredients 
ensuring that it exhibits both a finite-time singularity (and its associated scale invariance) and oscillatory behavior. 
The scale invariance is thus partially broken by the existence of dynamically generated length scales associated with 
the oscillations. We start from ([[]) and enrich it by the minimal ingredient to obtain what we believe is the simplest 
"normal form" of an oscillatory finite-time singularity. While the singularity emerges from the nonlinear growth law 
with positive feedback, the hierarchy of length scales results from a nonlinear negative feedback. The competition 
between the positive and negative nonlinear feedbacks create an approximate self-similar oscillatory structures, which 
can be understood from a spiral dynamics in phase space around a central unstable fixed point. Physically, the self- 
similar oscillations result from the dependence of the local frequency of the nonlinear oscillator on the amplitude. 
This will be shown in phase space to result from the special role played by the origin which is the unstable fixed point 
around which the spiral structures of trajectories are organized. 

This spiral structure of the dynamics around the central unstable fixed point bears a superficial resemblance to 
the the Shilnikov's mechanism for chaos [22]. However, both their dynamics and their behaviors are unrelated. 
Shilnikov's systems are characterized by trajectories in phase space spiraling towards the hyperbolic point along the 
stable manifold and then blowing-up exponentially along the unstable manifold of the hyperbolic point, until they are 
reinjected again along the stable manifold. In our system, trajectories in phase space spiral out slowly at first and then 
accelerate until a singular point in finite-time is reached due to a faster-than-exponential acceleration. Our system has 
thus a finite lifetime while Shilnikov's systems are globally statistically stationary. 

Our work is somewhat more related to that of several authors who emphasized the possible role of spiral structures 
in singular flows as a mechanism to promote the transfer of energy from large scales to small scales [Q, M|. Kiehn 
[34] has emphasized that vortex sheet evolution, governed by an integral form of the Biot-Savart law (known as the 
Rott-Birchoff equation) leads to the production of discontinuities in finite time. Asymptotic spiral type solutions 
in the vicinity of the singularity have been investigated both analytically and numerically (see Q] and references 
therein). Szydlowski et al. [54] have analyzed a nonlinear second-order ordinary differential equation, called the 
Kaldor-Kalecki business model in which capital stock changes are caused by past investment decisions. Their study 
emphasizes the negative feedback connected with the lag-delay effect and thus lacks the positive feedback trend effect 
discussed here. Canessa [^ [7|] has also a nonlinear second-order differential equation for the price but again the 
emphasis is on the nonlinear feedback rather than on the possibility of explosive phases coupled with the oscillatory 
behavior. 

Let us also mention another mechanism for log -periodicity: scale invariant equations which present an instability 
at finite wavevector decreasing with the field amplitude may generate naturally a discretely scale-invariant spectrum 
of internal scales [50]. 

We first motivate the normal form studied here for an oscillatory finite-time singularity by three physical examples, 
namely the time evolution of a stock market price described in section 2, the dynamics of human population described 
in section 3 and the coupled evolution of a damage variable and the average stress leading to material rupture given 
in section 4. We then present in section 5 an analysis of the effect of each of the two components (the nonlinear 
amplification and the nonlinear reversal term) of the dynamics taken separately. Section 6 describes in a rather 
heuristic way the fundamental characteristics of the overall dynamics obtained when combining both terms. Section 
7 provides a detailed dynamical system approach giving a complete characterization of the dynamics in phase space 
and precise predictions on the exponents of the scaling laws which are tested by numerical simulations. Section 8 
concludes. 
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2 Stock market price dynamics 



The importance of the interplay of two classes of investors, so-called fundamental value investors and technical 



analysists (or trend followers), has been stressed by several recent works Q40| , |lq ] to be essential in order to retrieve 
the important stylized facts of stock market price statistics. We build on this insight and construct a simple model of 
price dynamics, whose innovation is to put emphasis on the fundamental nonlinear behavior of both classes of agents. 

2.1 Nonlinear value and trend-following strategies 

The variation of price of an asset on the stock market is controlled by supply and demand, in other words by the net 
order size 17 through a market impact function [|i3|]. Assuming that the ratio p/p of the price p at which the orders 
are executed over the previous quoted price p is solely a function of 17 and using the condition that it is impossible 
to make profits by repeatedly trading through a close circuit (i.e. by buying and selling with final net position equal 



to zero), Farmer [15] has shown that the logarithm of the price is given by the following equation written in discrete 
form 

lnp(t + 1) - lnp(t) = ffi . (3) 

lj 

The so-called "market depth" L is the typical number of outstanding stocks traded per unit time and thus normalizes 
the impact of a given order size 17(i) on the log-price variations. The net order size 17 summed over all traders is 
changing as a function of time so as to reflect the information flow in the market and the evolution of the traders' 
opinions and moods. A zero net order size $7 = corresponds to exact balance between supply and demand. Various 
derivations have established a connection between the price variation or the variation of the logarithm of the price to 
factors that control the net order size itself [[TJ, |5[ |44fl. Two basic ingredients of 17(i) are thought to be important in 
determining the price dynamics: reversal to the fundamental value (17f un d(i)) an d trend following (17 trcn d(£))- Other 
factors, such as risk aversion, may also play an important role. 

We propose to describe the reversal to estimated fundamental value by the contribution 

^fund(i) = -c [mp(i) - hxpf] | hxp(t) - lnp/T" 1 , (4) 

to the order size, where p/ is the estimated fundamental value and n > is an exponent quantifying the nonlinear 
nature of reversion to p*. The strength of the reversion is measured by the coefficient c > 0, which reflects that 
the net order is negative (resp. positive) if the price is above (resp. below) pj. The nonlinear power law [lnp(i) — 
Inpf] | In p(t) — lnpy| n_1 of order n is chosen as the simplest function capturing the following effect. In principle, 
the fundamental value pf is determined by the discounted expected future dividends and is thus dependent upon 
the forecast of their growth rate and of the risk-less interest rate, both variables being very difficult to predict. The 
fundamental value is thus extremely difficult to quantify with high precision and is often estimated within relatively 
large bounds [|T], [R], 38, |||: all of the methods of determining intrinsic value rely on assumptions that can turn out 



to be far off the mark. For instance, several academic studies have disputed the premise that a portfolio of sound, 
cheaply bought stocks will, over time, outperform a portfolio selected by any other method (see for instance [[37|]). As 
a consequence, a trader trying to track fundamental value has no incentive to react when she feels that the deviation 
is small since this deviation is more or less within the noise. Only when the departure of price from fundamental 
value becomes relatively large will the trader act. The relationship ([|) with an exponent n > 1 precisely accounts 
for this effect: when n is significantly larger than 1, \x\ n remains small for |x| < 1 and shoots up rapidly only 
when it becomes larger than 1, mimicking a smoothed threshold behavior. The nonlinear dependence of 17f un d(t) on 
ln[p(t)/pj] = \np(t) — Inpf shown in (Q) is the first novel element of our model. Usually, modelers reduce this term 
to the linear case n = 1 while, as we shall show, generalizing to larger values n > 1 will be a crucial feature of the 
price dynamics. In economic language, the exponent n = dm 17f un d/<im (ln\p(t)/pf]) is called the "elasticity" or 
"sensitivity" of the order size 17f un d with respect to the (normalized) lOg-price ln\p(t)/pf]. 

A related "sensitivity", that of the money demand to interest rate, has has been recently documented to be larger 
than 1, similarly to our proposal of taking n > 1 in (H). Using a survey of roughly 2,700 households, Mulligan 
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and Sala-i-Martin [43] estimated the interest elasticity of money demand (the sensitivity or log-derivative of money 
demand to interest rate) to be very small at low interest rates. This is due to the fact that few people decide to 
invest in interest-producing assets when rates are low, due to "shopping" costs. In contrast, for large interest rates 
or for those who own a significant bank account, the interest elasticity of money demand is significant. This is a 
clear-cut example of a threshold-like behavior characterized by a strong nonlinear response. This can be captured by 
e = dhxM/dhxr = (r/r in fl) n with n > 1 such that the elasticity e of money demand M is negligible when the 
interest r is not significantly larger than the inflation rate r in Q and becomes large otherwise. 

Trend following (in various elaborated forms) was (and probably is still) one of the major strategy used by so- 
called technical analysts (see ^ for a review and references therein). More generally, it results naturally when 
investment strategies are positively related to past price moves. Trend following can be captured by the following 
expression of the order size 

^trend(t) = ai[lnp(t) - lnp(i - 1)] + a 2 [lnp(t) - ]np(t - 1)]| lnp(t) - \np(t - l)^ 1 . (5) 

This expression corresponds to driving the price up if the preceding move was up (ai > and a 2 > 0). The linear 
case {a\ > 0, a 2 = 0) is usually chosen by modelers. Here, we generalize this model by adding the contribution 
proportional to a 2 > from considerations similar to those leading to the nonlinear expression (Q) for the reversal 
term with an exponent n > 1. We argue that the dependence of the order size at time t resulting from trend-following 
strategies is a nonlinear function with exponent m > 1 of the price change at previous time steps. Indeed, a small 
price change from time t — 1 to time t may not be perceived as a significant and strong market signal. Since many of 
the investment strategies are nonlinear, it is natural to consider an average trend-following order size which increases 
in an accelerated manner as the price change increases in amplitude. Usually, trend-followers increase the size of their 
order faster than just proportionally to the last trend. This is reminiscent of the argument ^ that traders's psychology 
is sensitive to a change of trend (acceleration or deceleration) and not simply to the trend (velocity). The fact that 
trend-following strategies have an impact on price proportional to the price change over the previous period raised to 
the power m > 1 means that trend-following strategies are not linear when averaged over all of them: they tend to 
under-react for small price changes and over-react for large ones. The second term with coefficient a 2 captures this 
phenomenology. 

2.2 Nonlinear dynamical equation for stock market prices 

Introducing the notation 

x(t) = hx\p{t)/p f ] , (6) 
and the time scale St corresponding to one time step, and putting all the contributions (|j) and (|5|) into (|3|), with 

D,(t) = fifund(i) + ^trend(i), We get 

x(t + St)-x(t) = j (ai [x(t) - x(t - St)] + a 2 [x(t) - x(t - 6t)]\x(t) - x(t - St^ 171 " 1 - c x(t)\x(t)\ n -^ . (7) 
Expanding ([7J) as a Taylor series in powers of St, we get 



4 ,* +fa ^* | * r . 1 _- l( , )|i „ )r , H cW| . 



where O [(5i) 3 ] represents a term of the order of (St) 3 . Note the existence of the second order derivative, which results 
from the fact that the price variation from present to tomorrow is based on analysis of price change between yesterday 
and present. Hence the existence of the three time lags leading to inertia. A special case of expression (|7|) with a 
linear trend-following term (a 2 = 0) and a linear reversal term (n = 1) has been studied in [Q, [15|], with the addition 
of a risk-aversion term and a noise term to account for all the other effects not accounted for by the two terms (Q) and 
(Q). We shall neglect risk-aversion as well as any other term and focus only on the reversal and trend-following terms 
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previously discussed to explore the resulting price behaviors. Grassia has also studied a similar linear second-order 
differential equation derived from market delay, positive feedback and including a mechanism for quenching runaway 



markets [19]. Thurner [55] considers a three-dimensional system of three ordinary differential equations coupling 
price, "friction" and a state variable controlling friction, which can be mapped onto a third-order ordinary differential 
equation. The nonlinearity is on the friction term and not on the trend term which is again assumed linear. 

Expression (0) is inspired by the continuous mean-field limit of the model of Pandey and Stauffer [Q], defined by 



starting from the percolation model of market price dynamics [N4L 12, 53] and developed to account for the dynamics 



of the Nikkei and Russian market recessions [26, ^7|]. The generalization assumes that trend-following and reversal 
to fundamental values are two forces that influence the probability that a trader buys or sells the market. In addition, 



Pandey and Stauffer [44] consider as we do here that the dependence of the probability to enter the market is a 
nonlinear function with exponent n > 1 of the deviation between market price and fundamental price. However, they 
do not consider the possibility that m > 1 and stick to the linear trend-following case. We shall see that the analytical 
control offered by our continuous formulation allows us to get a clear understanding of the different dynamical phases. 

Among the four terms of equation (|J), the first term of the right-hand-side of (||) is the least interesting. For 
at < L, it corresponds to a damping term which becomes negligible compared to the second term in the terminal 
phase of the growth close to the singularity when \ dx/dt\ becomes very large. For a\ > L, it corresponds to a negative 
viscosity but the instability it provides is again subdominant for m > 1. The main ingredients here are the interplay 
between the inertia provided by the second derivative in the left-hand-side, the destabilizing nonlinear trend-following 
term with coefficient a2 > and the nonlinear reversal term. In order to simplify the notation and to simplify the 
analysis of the different regimes, we shall neglect the first term of the right-hand-side of ((8)), which amounts to take 
the special value a\ = L. In a field theoretical sense, our theory is tuned right at the "critical point" with a vanishing 
"mass" term. 

Equation (||) can be viewed in two ways. It can be seen as a convenient short-hand notation for the intrinsically 
discrete equation (0), keeping the time step St small but finite. In this interpretation, we pose 

a = a 2 (St) m - 2 /L, (9) 
7 = c/L(St) 2 , (10) 

which depend explicitely on St, to get 

A second interpretation is to genuinely take the continuous limit St —>■ with the constraints a^/L ~ (St) 2 ~ m and 
c/L ~ (St) 2 . This allow us to define the now <5t-independent coefficients a and 7 according to @ and (p"l|) and obtain 
the truly continuous equation (|ll|). This equation can also be written as 

~& = 2/2 ' (12) 

= «2/2|2/2| m - 1 -7yi|yir 1 - (13) 

at 

This is the system we are going to study for m > 1 and n > 1. For further discussions, we call the term proportional 
to a (resp. 7) the trend or positive feedback term (resp. the reversal term). The richness of behaviors documented 
below results from the competition between these two terms. 

In defining the generalized dynamics ( |I2| , |T"3| ) for the market price, we aim at a fundamental dynamical understand- 
ing of the observed interplay between accelerating growth and accelerating (log-periodic) oscillations, that have been 
documented in speculative bubbles preceding large crashes [31, ^7|, [2~8| ]. 



We shall show below that the origin (y\ = 0, 1/2 = 0) plays a special role as the unstable fixed point around which 
spiral structures of trajectories are organized in phase space (2/1,2/2)- It is particularly interesting that this point plays 
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a special role since y\ = means that the observed price is equal to the fundamental price. If, in addition, yi = 0, 
there is no trend, i.e., the market "does not know" which direction to take. The fact that this is the point of instability 
around which the price trajectories organize themselves provides a fundamental understanding of the cause of the 
complexity of market price time series based on the instability of the fundamental price "equilibrium". 



3 Population dynamics 

As a standard model of population growth, Malthus' model assumes that the size of a population increases by a fixed 
growth rate a independently of the size of the population and thus gives an exponential growth: 

(i4) 

The logistic equation attempts to correct for the resulting unbounded exponential growth by assuming a finite carrying 
capacity K (t) such that the population instead evolves according to 

± = a p(t)[K(t)-p(t)] , (15) 
where <tq controls the amplitude of the nonlinear saturation term. Applying this model to the human population 



on earth, Cohen and others (see [13] and references therein) have put forward idealized models taking into account 
interaction between the human population p(t) and the corresponding carrying capacity K(t) by assuming that K(t) 
increases with p(t) due to technological progress such as the use of tools and fire, the development of agriculture, 
the use of fossil fuels, fertilizers etc. as well an expansion into new habitats and the removal of limiting factors by 
the development of vaccines, pesticides, antibiotics, etc. If K(t) grows faster than p(t), then p(t) explodes to infinity 
after a finite time creating a singularity [[3C|]. In this case, the limiting factor — p(t) can be dropped out and, assuming 
a simple power law relationship K(t) oc [p(t)] <5 with 6 > 1, (15) can be written as ( |l"4| ) with an accelerating growth 
rate a replacing a$: 

a = a [p(t)} 5 . (16) 

The generic consequence of a power law acceleration in the growth rate is the appearance of singularities in finite 
time: 

p{t) = p(0) — , for t close to t c , (17) 



t c 

where t c is determined by the constant of integration, i.e., the initial condition p(0) as t c = [p{0)]~ s /Sao. Equation 
( |l6| ) is said to have a "spontaneous" or "movable" singularity at the critical time t c [Q], 



Note that, using (|17|), ( |16| ) can be written 

da 



aB >, (18) 



showing that the finite-time singularity of the population p(t) is the result of the finite-time singularity of its growth 
rate, resulting from the quadratic growth equation (|I8|). 
We now generalize (18) as 

da 



aa\a\ m ^ - 1 \n{p/K OQ )\\n{p/K 00 )\ n - 1 (19) 
dt 

for the following reasons. Apart from the absolute value, the first term in the r.h.s. of ( |T9| ) is the same as (18) with 
m = 2. In addition, we allow the instantaneous growth rate a to be negative and thus its growth has to be signed. 
The novel second term in the r.h.s. of ( ff9| ) takes into account a saturation or restoring effect such that by itself this 
term attracts the population p(t) to an asymptotic constant carrying capidity K^. Using the logarithm of the ratio 
pjKoa is the natural choice for the dynamics of a growth rate since hxip/K^) is nothing but the effective cumulative 
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growth rate. For n = 1, — 7 \n.{p/ ln(p/ K^Y" 1 = — r yln{p/K 00 ) corresponds to a linear (in \n(p/ K^)) 
restoring term. A choice n > 1 captures the following effect: the restoring term is very weak when p departs weakly 
from K and then becomes rather suddenly stronger when this deviation increases. This nonlinear feedback effect is 
intended to capture the many nonlinear (often quasi-threshold) feedback mechanisms acting on population dynamics. 
In the limit n — > +00, the reversal term acts as a threshold. Note that the absolute values can be removed when the 
exponents m and n are odd. 

Expression (19) generalizes ( |l"5| ) by putting together a faster-than-exponential growth and an attraction to finite 
value. In contrast, (|15|) puts together an exponential growth and an attraction to a finite value. 
Let us introduce the change of variable 

yt = Hp/K^), (20) 

2/2 = 0" • (21) 

The equations (14,19) then retrieve (|l2|) and (13). For further discussions, we shall refer to the term proportional to a 
(resp. (3) as the positive feedback or acceleration term (resp. the reversal term). 

In defining the generalized dynamics (12) and (13) with ( [20] ) for the population, we aim at a fundamental dynam- 
ical understanding of the observed interplay between accelerating growth and accelerating (log-periodic) oscillations, 
that have been documented in [33 , 5^, I 



,30] 



4 Rupture of materials with competing damage and healing 



Consider the problem of so-called creep or damage rupture [ ]35[ ] in which a rod is subjected to uniaxial tension by 
a constant applied axial force P. The intact cross section A(t) of the rod is assumed to be a function of time. The 
physical picture is to envision myriads of microcracks damaging progressively the rod and decreasing its effective 
intact cross section that can sustain stress. The problem is simplified by assuming that A(t) is independent of the 
axial coordinate, which eliminates necking as a possible mode of failure. The considered viscous deformation is 
assumed to be isochoric, i.e., the rod volume remains constant during the process. This provides a geometric relation 
between the rod cross-sectional area and length AqLq = A(t)L(t) = constant, which holds for all times. 
The rate of creep strain e c can be defined as a function of geometry as 

de c _ 1 dL _ 1 dA 
~dt ~ L~dJ ~ A~dJ 

showing that 



(22) 



e^Hln^-lnf), (23) 

where Lq = L(to) and Aq = A(to) correspond to the underformed state e c (io) = at time to- 

The rate of change of the creep strain e c (i) is assumed to follow the rheological Norton's law, i.e., 

^ = Ca» , with ii > , (24) 
dt 

where the stress 

a = Pj A (25) 



is the ratio of the applied force over the cross section of the rod. Eliminating de c /dt between (22) and ([P|) and 

using @ leads to A^dA/dt = -CP^, i.e., A(t) = A(0) (^i^) 1 ^', where the critical failure time is given by 
t c = [A(0) /PY 'I il^C). The rod cross section thus vanishes in a finite time t c and as a consequence the stress diverges 
as the time t goes to the critical time t c as 
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Physically, the constant force is applied to a thinner cross section, thus enhancing the stress, which in turn accelerate 
the creep strain rate, which translates into an acceleration of the decrease of the rod cross section and so on. In other 
words, the finite-time singularity results from the positive feedback of the increasing stress on the thinner cross section 
and vice-versa. This finite-time singularity for the stress can be reformulated as a self-contained equation expressed 
only in terms of the stress: 

Tt = Ca " +1 ■ (27) 
Let us now generalize this model by allowing not only creep deformations leading to damage but also recovery or 



aa\a 



m—l 



7e c |e, 



m-l 



(28) 



healing as well as a strain-dependent loading. We thus propose to modify the expression (27) into 

da 
~dt 

The first term in the right-hand-side of (^) is similar to (|27|) by redefining fj, + 1 as m, and captures the accelerated 
growth of the stress leading to a finite-time singularity. It embodies the positive geometrical feedback of a reduced 
intact area on the effective stress applied to whole system. The addition of the second term in the right-hand-side of 



(28) implies a modification of Norton's law which is no more specified by the exponent fiorm and introduces the 
novel physical ingredient that damage can be reversible. For convenience, we choose a specific power law dependence 
— 7e c |e c | n_1 to capture the healing mechanism. This term alone tends to decrease the effective stress and describes a 
recovery of the material since a reduction of the effective stress is associated with an increase of carrying area of the 
intact material. Alternatively, we can interpret (|28|) as defining the loading, which becomes strain-dependent: a larger 
strain implies less room for additional stress increase, as for instance occurs in mechanical apparatus in the laboratory 
which are often limited to small deformations and relax the applied stress beyond a given strain. The mechanism is 
also attractive for describing the tectonic loading of faults which is occurring with mixed stress and strain rates, rather 
than a pure imposed stress or strain rate. 



Bringing the system out of equilibrium and then releasing it, the equation (28) describes how the system can 
either recover an equilibrium or rupture in finite-time due to accumulating creep and damage in its dynamical attempt 
to come back to equilibrium. The novel second term in the r.h.s. of (^8|) takes into account a healing process or 
work-hardening, such that large creep deformations hinder and may even reverse the stress increase. By itself, this 
term attracts the cross section A(t) back to the equilibrium value Aq. Since the cross-sectional area A(t) can be 
alternatively interpreted as the surface of intact material able to carry the stress, healing increases the area of intact 
material and thus decreases the effective stress. 

We close the model by assuming again Norton's law but with an exponent fjf different from ji: 



^ = Co* 
dt 



(29) 



Incorporating the constant C in a redefinition of time Ct 
and j/C — > 7) and posing 



Vi 

y-2 



t (with suitable redefinitions of the coefficients a/C —> a 



In 



An 



(30) 
(3D 



we retrieve the dynamical system (12) and (13) for the special choice // = 1, which we shall restrict to in the sequel. 
We are going to study this system in the regime where m > 1 and n > 1. The first condition m > 1 ensures the 
existence of a finite-time singularity describing a positive feedback between the stress increase and the cross section 
decrease. The second condition n > 1 ensures that the healing process is only active for large deformations: the larger 
n is, the more threshold-like is this effect with respect to the amplitude of the creep strain. 

In defining the generalized dynamics (12(13) with (30) for the rupture dynamics, we aim at a fundamental dynam- 
ical understanding of the observed interplay between accelerating growth and accelerating (log-periodic) oscillations, 
that have been documented in time-to-failure analysis of material rupture [^, 57, 2[ 47, 17, 21, 25, [29| ]. 
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5 Individual components of the dynamics 



The system ( |T2[ , |l3| ) can also be written as (11), which we rewrite here for convenience, with uniform notation: 

d 2 yi 



dt 2 



i . dy x dy x x 



(32) 



This autonomous expression has the following interpretation. The left-hand side (l.h.s.) is the inertia for the variable 
y\. The right-hand side (r.h.s.) has two elements. The first element in the r.h.s. describes a restoring force for 7 > 
which will be the case studied here. Coupled with inertia, this leads generally to an oscillatory behavior. The second 
element is a positive trend which can lead to singular behavior for a > and m > 1. 

Our strategy is first to study these two components coupled to the inertia as separate sub-dynamical systems. The 
interplay between these two components will then help deepen our understanding of the overall dynamics. 

Each sub-dynamical system can be reduced to a normalized model where one single (resp. pair of) orbit(s) 
give(s) a template dynamics of the entire sub-dynamical system. Because the system is autonomous, we describe 
two complementary approaches to describe its dynamics: (1) phase portrait using orbits graphed in the phase space 
y = (yi, 2/2); (2) time evolution of trajectories y(t; yo, to) with initial condition yo at time to- 

5.1 The restoring term: oscillations 

We first consider the restoring term and examine the nature of the corresponding oscillations. Motivated by sections 
||, H and ^| describing the application to various physical processes, we focus on the relevant case 7 > and n > 1. 



5.1.1 Model 

Keeping only the first term in the r.h.s. of ([32|) gives: 

dt 2 



-7 2/1I2/1 



n-1 



This system can be written as a one-degree-of-freedom Hamiltonian system 

d / 2/1 \ _ / 3^ TJU. \ _ / 2/2 



dt \ 2/2 



dyi 



"72/1 1 2/1 



ln-1 



where 



n + l 



o n + l 1 r, 

(2/?) 2 +^y 2 2 - 



(33) 



(34) 



(35) 



An orbit of ( |34| ) in the y phase space for fixed (n, 7) can be given as a graph of constant H (y; n, 7). In other words, 
a trajectory y(t; yo, to) going through yo follows contours of constant H(yo; n, 7) in time. 



5.1.2 Template dynamics in the normalized model 

We now describe the template dynamics using a normalized model. We define the following normalized variables 
denoted by hat {•}: 



1 2/1 

2/2 



( 7"+i m \ 



2/1 

2/2 



(36) 
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Then all H(y; n, 7)-contours in the original y phase space collapse onto a single closed curve in the normalized y 
phase space. The case H = at y = 0, which is an elliptic fixed point, is excluded from this analysis. The closed 
curve corresponds to the normalized orbit defined by: 

1 , .o, n + 1 1 



1 



n + 

Therefore, the range of oscillation for the variable y\ and its velocity y 2 = J^yi is: 

yx € [-(n+l)^,(n + l)^] , y 2 e[-y/2,V2], (38) 

as shown in the left panels of Figure |[ The original dynamics in y can be recovered using (36). 

The closed orbit in phase space represents an oscillatory dynamics for any fixed n > 0. For n = 1, the curve given 
by (37) is a perfect circle with radius y/2 corresponding to a linear harmonic oscillator (Figure [[]:). For n / 1, the 
closed curve is a nonlinear generalization of the circle. As n increases above 1, (yf ) 2 becomes small for \yi\ < 1. 
From the condition (37), it follows that \y 2 \ remains almost constant near its maximum \/2 for \y\ \ < 1. Furthermore 
from (38), the range of |yi| decreases monotonically to 1 as n increases, while the range of y 2 remains unchanged for 
any n. On the whole, the geometry of the closed orbit becomes closer and closer to a square as n increases, as shown 
in Figures [T]e and [j]g. 

The normalized dynamical model can be obtained by substituting ([3f]) into (34). In addition, it can be decoupled 
into two independent first-order nonlinear ordinary differential equations (ODEs) using the condition (|3 

m \_( sign[y 2 ] V2 [1 - Ax{Vx)^\ u 




-yx\yi\ n x J \-siga[yx] (n + l)^TT [1 _ U$\ 



) ■ (39) 



n+1 v 

h 

Therefore, the normalized trajectory goes around the origin in a clockwise direction along the closed curve given by 
(|37|). 



Given an initial condition yo satisfying (pl\), the two ODEs in (g9j) can be solved separately and provide the 
trajectory y(t;yo,£o)- However, once yi(i;yo,^o) is solved using the top ODE, 2/2(^5 yci,io) can be algebraically 
computed using ([37|) without solving the bottom ODE for y 2 ; and vice versa. The period of oscillation can be obtained 
by integrating one of the two ODEs over the range defined by (38). 

For some special re's, the normalized model has explicit analytical solutions. For n = 1, corresponding to a 
harmonic oscillator, the solution is: 



y(*;yo,io) = V2(sin(t-to + ^),cos(t-t + ^)) , (40) 

where the initial condition yo = \/2(sin 6, cos 9) satisfies ( |37[ ) for any real number 9. Its period of oscillation is 2ir. 
For n = 3, the solution for yi(i; yo, to) is Jl8|]: 



F ( arccos -±=) = y w (t - t ) , (41) 

where yio = y\(t = £0) and fotf = to) = and F (u, k)) is the elliptic integral of the first kind defined by 

f u fill 

J o v 1 — fe sin v 

The evolution of a trajectory with initial condition yo = (0, v2) at to = is shown in the right panels of 
Figure [I] as a function of i. For n = 1, the normalized variable y\ is preceded by its velocity §2 by a ^-phase 
shift. As n increases, the amplitude of the normalized velocity I2/2I becomes nearly constant about v2 for \yx\ < 1 
[see Figures |I]e and [l^ as well as the discussion concerning the geometry of the normalized orbit given by (37)]. 
Accordingly, y\ evolves almost linearly in time with a constant velocity y 2 = ±y/2 for |yi| < 1. Thus, time series 
of yi and y 2 respectively have saw-teeth and step-function shapes. The period of the oscillation becomes shorter as 
n increases, because the velocity amplitude \y 2 \ remains closer and closer to its maximum as n increases during each 
oscillation cycle. 
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5.1.3 Global dynamics 



Having understood the template dynamics of the normalized model, we describe the global dynamics of the oscillatory 
element in the original model. Two parameters, n and 7, are used to define the one-degree-of-freedom Hamiltonian 
system given by From the normalization defined by (36) and the resulting system given by (|39|), we see that the 
exponent n is the controlling parameter. The coefficient 7 contributes only for the scaling of (yi,t). 

For fixed (n, 7), the phase portrait in y consists of a family of closed orbits parameterized by H. Each orbit has a 
unique H and oscillates in the clockwise direction around the origin (left panels of Figure ^|). Given the properties of 
the normalized orbit in Section 5.1.2 , three main properties of individual orbit as function of H follow. 

First, H measures the amplitude of oscillation. Each orbit in y ranges over: 



2/1 6 



(n + l)H 

7 



(n + l)H 



1 



y 2 G [-(2H)2,(2H) 



(43) 



from ( p6[ ) and (|38|). The left panels of Figure g show the orbits in the y phase space as curves of constant H. The 

higher H is, the larger is the amplitude of the oscillations. 

Second, H determines the geometry of the orbits, except for n = 1 where all orbits are ellipses of the same aspect 
1 _ 

ratio 7 ™+i (Figure Nc). For n / 1, we describe the geometry in terms of the deformation from the normalized curve 

r- — 1 . . _L_ 1 

( p7| ) by comparing the two normalization coefficients for y\ and y 2 in ([4 3D, i.e., H n + l and H 2 , respectively. These two 

coefficients also govern the range of the oscillation (36). For n < 1, the two coefficients have the following relation: 

1 1 _j_ 1 

iJn+i < H2 < 1 for H < 1, and ff «+i > H2 > 1 for H > 1. Hence, an orbit with small amplitude (H < 1) has 

a geometrical shape stretched along the vertical direction (Figure ||a), as deduced from the normalized closed curve 

(Figure [TJa). This is because y\ is reduced more than y 2 . Similarly, an orbit with a large amplitude (H > 1) has a 

_j_ 1 

geometrical shape stretched horizontally. For n > 1, the relation is reversed: 1 > H n + l > H2 for H < 1, and 

1 1 

1 < H "+ 1 < H 2 for H > 1. Accordingly, orbits with small (H < 1) or large (H > 1) amplitudes respectively have 
horizontally or vertically stretched geometries compared with the normalized closed curve (Figure |2|c,d). 

Third, H controls the speed of the oscillation, except for the harmonic oscillator case n = 1 which has a constant 
period 2^/^/7. For n/1, the period of the oscillations is obtained using (p6|), (38) and (39): 



T(H;n,j) 
where C(n) is a positive number given by: 

C(n) 



C(n) 7" 



v/2 



1— n 
,H~2(n + l) 



(n + i; 



n + 1 JO 



dy 2 



(1 - 



(44) 



(45) 



n+l 



Differentiating the expression of the period T(H; n, 7) given by (44) with respect to H gives: 



^T(E T;„, 7 ) 



1 



3n + l 



2(n + i; 



C(n) 7™+! ^+T) 



(46) 



Therefore, -^jT(H;n,j) can be positive or negative depending on n. For n < 1 where > 0, the period 

increases monotonically from to 00 as H increases. In contrast, for n > 1, the period decreases monotonically from 
00 to as H increases. The right panels of Figure || show the period of oscillation on the abscissa as function of the 
maximum amplitude reached by y 2 equal to y/2H, given as a measure of the oscillation amplitude. 



5.2 The trend term: singular behavior 

We now consider the trend term and examine the nature of the singularity which manifests itself in "finite-time" in 
the behavior of the velocity y 2 and of the variable y\ . Motivated by the applications to concrete physical processes 
discussed in sections |2[ ^]and |], we focus on the case a > and m > 1. 
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5.2.1 Model 



Keeping only the second term in the r.h.s. of (32) and re-writing it as a system of two-dimensional ODEs for the 
variable y\ and its velocity y 2 give: 

d I yi \ I m \ 



dt \ y 2 



\m— 1 



In this system, 



G(y;m,a) = yi 



*(2-m) 



2/2I2/2I 

sign [y 2 ] log |y 2 | 



for m ^ 2 
for m = 2 



(47) 



(48) 



is an invariant, i.e., ^G(y; m, a) = 0. This can be easily verified by differentiating ( J48| ) with respect to t and then 
substituting in (p7]). Therefore, like H(y; n, 7) in the one-degree-of-freedom Hamiltonian system discussed in the 
previous section for the restoring term, a curve of constant G(y; m, a) in the y phase space corresponds to an orbit 
governed by @. A trajectory y(t; y ,t ) going through y satisfies G(y(t; y , t ); m, a) = G(y ; m, a) for any t. 



5.2.2 Template dynamics in the normalized model 

To describe the template dynamics of the trend term, we reduce the system ( p7| ) to a normalized model using the 
invariant G (48). We define the following normalized variables denoted by {•}: 



( VI 

2/2 



t 



{ a 



III 

112 



+ G 



(49) 



a 



Substituting ( |49[ ) into ( |47| ) gives the normalized dynamical model: 




2/2 



Therefore, the normalized velocity y 2 undergoes an irreversible amplification if it starts from |y 2 
the variable y\ along. The dividing point y 2 = is a fixed point. Furthermore, substituting (49) into 
invariant condition for the normalized orbit: 



(50) 

/ and it pulls 
gives an 



2/1 = K2/2) 



(51) 



where 



(2-m) 

sign[y 2 ] 



„~ |1— m 
2/2 1 2/2 1 

log 1 2/2 1 



for m 7^ 2 
for m = 2 



The left panels of Figure || show graphs of (^Tj) for m = 1.5, 2 and 2.5. Each graph consists of a pair of 
rotationally-symmetric orbits for y 2 > and y 2 < with the symmetry given by y — > — y. Through the normalization 
d49|), all orbits in the original y phase space with y 2 > collapse onto the single normalized orbit with y 2 > in 
y phase space. Similarly, all orbits in the y phase space with y 2 < collapse onto the single normalized orbit with 
y 2 < in the y phase space. Using (|50|), the slope of the graph is: 



dm = dy 2 ,dyi 
dyi dt dt 



\m\ m - 1 



(52) 



Integrating this equation retrieves the normalized invariant condition given by d5~T|). For fixed m, the slope increases 
from to 00 as \y 2 \ increases. The higher m is, the faster the slope increases. 
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The geometry of (the pair of) symmetric normalized orbit(s) depends on m significantly (Figure |J). This is due 
to a qualitatively different behavior of y\ = (7(2/2) as 2/2 approaches towards either end of the normalized orbit. For 

2/2 > 0, 

for 1 < m < 2 

2/1 = 9(2/2) G <( (-00,00) form = 2 , y 2 e (0, 00) . (53) 

for m > 2 

In other words, the normalized model undergoes a bifurcation at m = 2. Precisely for m = 2, 9(2/2) spans the whole 
interval (— 00, 00) and thus interpolates between the cases for m < 2 and m > 2 for which only one-half of the 
interval is covered by the dynamics. 

Accordingly, the dynamical behavior of the normalized trajectory changes qualitatively when 1 < m < 2, m = 2 
or m > 2. Any initial condition yo = (2/1,0,2/2,0) °f the normalized trajectory y(t;y ,to) must satisfy (^T|), i.e., 
2/1,0 = 9(2/2,0)- Using this initial condition, ( |50| ) can be solved analytically: 

y\{h yo> *o) = 9(2/2,0) + sign[2/ 2 ,o] 

( m _ 1)— _J_ [(t c (y 2j0 ) - t)— - (t c (y 2 ,o) - t )~] for m ^ 2 
log f r#4s) for m = 2 ' 

1 _ ^ 1 

2/2(*;yo,*o) = sign[y 2 ,o] (m - 1) {t c (y 2 ,o) ~ t) . (54) 
The solution is valid only for a semi-infinite time interval t G (—00, t c (2/2,o)) up to the normalized "critical time:" 

t c (m,o) = k + —^—Imfi] 1 -™, (55) 

m — 1 

for any arbitrary to. The center panels of Figure ^] shows t c on the abscissa as a function of 2/2,0 on the ordinate with 
to = 0. The solution becomes singular as 2/2 approaches ±00, and no solution exists for t > i c (?/2,o)- 

This is the "finite-time" singular behavior for the velocity amplitude \y 2 \, which occurs for any m > 1. It is driven 
by the nonlinear positive feedback of the trend term producing a faster than exponential growth rate, leading to a 
infinite growth of | y 2 \ in finite time. For a fixed m, the initial normalized velocity 2/2,0 solely determines the behavior 
of the trajectory given by @ and @, as a consequence of the fact that the dynamics (|50|) is solely determined by 



y 2 . The evolution of y2(i; yo, to) in time with a pair of initial conditions yo = (0, ±0.6) at to = is shown in the 
right panels of Figure [| 

We now examine separately the behavior for forward (t > to) and backward (t < to) time intervals. Using the 
parity symmetry, we focus on the dynamics described by the normalized orbit with 2/2,0 > 0. Similar results hold for 
2/2,0 < using y -> -y. 

Over a forward finite time interval t G [to, t c (y2,o)) U P to ic(y2,o)> the normalized trajectory with initial condition 
2/2,0 > ranges over: 

I [9(2/20), 00) for 1 < m < 2 _ r _ . 
^(S,o)) form>2" • WSfocoo), ( 56 ) 

as shown in the left and right panels of Figure ||[ During this finite time interval up to to(2/2,o)» 2/2 grows from 2/2,0 
to 00. For 1 < m < 2, y\ also blows up to infinity, dragged along by y 2 - On the contrary for m > 2, 2/1 culminates 
at a finite value, given by 9(2/2,0)- This finiteness of 2/1 occurs only for m > 2 for which the relative growth rate of 



2/2 compared to 2/1 also becomes singular, as given by the slope of the graph ^jfe = ^|r/-Jr i n (52). The finite final 



value 9(2/2,0) °f 2/1 can t> e observed in Figure gg using the graph of the normalized orbit ( pT| ) by replacing (2/1,2/2) by 
(9(2/2,0), 2/2,o)- 
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In contrast, over a backward semi-infinite time interval t < to, the normalized trajectory ranges over: 

f (0,5(2/20)] for 1 < m < 2 „ fn _ , 
2/iG<^ v 'l)r ,u , . , 2/2 € (0, 2/2,0 ■ (57) 
[ (-00,5(2/2,0)] for m > 2 

The normalized velocity 2/2 hence shrinks to 0, resulting in a finite increment equal 7/2,0 over the whole time interval. 
For 1 < m < 2, the increment 5(2/2,0) in 2/1 is a l so finite. On the contrary for m > 2, the increment in y\ is infinite 
because the relative growth rate given by slope ^ diverges for 2/2 <C 1. 

5.2.3 Global dynamics 

We now derive the global dynamics from the template dynamics of the normalized model constructed with the trend 
term. Two parameters, m and a, are used to define the two-dimensional representation of the trend term given by d47| ) 



and (48). For fixed (m, a), the model has the parity symmetry: y — > — y. We refer to the pair of symmetric orbits 



corresponding to the graph of G(y; m, a) = as the "reference orbits". From the normalization Q, we see that the 
dynamics along the pair of reference orbits is related to the dynamics along the normalized orbits (|5l1)-(57) through 
scaling of (2/1, i) by aT x . We also see that all orbits in the y phase space collapse onto the corresponding reference 
orbit through linear translations in the y\ phase space given by the distance G(y;m,a). Therefore, the dynamics 
along any orbit is exactly the same as the one along the corresponding reference orbit, except for a translation in the 
yi phase space. 

Accordingly, the phase portrait for any m consists in a pair of symmetric families of open orbits. Each family is 
parameterized by G (Figure |]) where the reference orbits are labeled by "0" [see also the left panels of Figure || for 
comparison with the normalized orbits]. On y 2 = 0, the dynamics is at rest. Therefore, the y phase space can be 
divided into dynamically distinct regions as follows. 

Definition 5.2.1 Singular basins and and boundary 6 a i ng 
We define a pair of singular basins: 

^stng = {y| 2/2 > 0, where —y 2 > 0} 

^smg = {Yi 2/2 < 0, where ^y 2 < 0} . (58) 
The boundary which separate the two basins is defined by: 

d 
dt 

Any point yo in the phase space belongs to either one of -B^ ng , B~ ing or 6 s i ng , as shown in figure |5[ 
For any m > 1, the individual trajectory satisfies the following. 

Corollary 5.2.2 

Consider a trajectory y(t; yo, io) with initial condition yo = (2/1,0, 2/2,0) at ti me *o- If it starts in B^ ng , then it will 
remain within the basin without never leaving it. It will reach the corresponding finite-time singularity, i.e., 

if yo e £ s 1ng> then y(*; yo^o) e £ s t„ g > 

for t e (-00, t c (y 2 ,o)) with y 2 {-oo; y , to) = ±0 and 2/2(^(2/2,0); yo, *o) = ±00. (60) 

The critical time t c (2/2,o) = a_1 ^c(2/2,o) i s a function of the initial velocity y 2 ,o only. In the case where y(t; yo, *o) 
starts on 6 s i ng , it will remain on the boundary basin without never leaving it for a bi-infinite time interval, i.e., 

if yo easing, then y(t; y , t ) = y G 6 si ng, for t G (—00, —00) . 



^sing = {y|2/2 = 0, where — y 2 = 0} . (59) 
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Definition 5.2.3 Source strips 5 s | ng and S si 



'sing 



Right next to the boundary 6 s i ng in each basin -B^ ng , we define a thin vertical strip of constant width: 



S: 



± 

sing 



{y g b? 



sing 



1 7/2 1 < 1, where \fy 2 



1 2/2 



m+l 



«1} 



(61) 



where y(i; yo, to) moves away from 6 s j ng extremely slowly. Once it leaves , it then reaches very quickly the finite- 
time singularity. For any yo, integrating with backward time, any y(t; yo, to) Wli l eventually enter S^ ng . Therefore, 
5^ ng can be considered as source regions of the finite-time singularity. 



However, the qualitative behavior of each individual orbit in each singular basin depends on the specific value of 



the exponent m, as discussed in Section p.2.2 



6 Overall dynamics: Fundamental characteristics 
6.1 Normalized model 

We now consider the overall dynamics obtained by combining the restoring and trend terms which have been analyzed 
separately in Section |5[ We use the following normalized variables: 



(62) 



to minimize the number of parameters by removing the coefficient a of the positive feedback trend term. For simplicity 
in the notations, we drop the bar {•} from here on. Then, the overall dynamical systems is written as: 



1 yi 


\ 


/ 


a 


yi 


\ 


2/2 








y-2 




t 






a 


t 




V 7 


/ 


V 


a -("+l) 


7 


/ 




y2 

2/2 osc + 2/2 sing 



where 



2/2 osc 
2/2 sing 



-72/1 1 2/1 

V2\V2\ m 



n-1 



(63) 



(64) 
(65) 



are the oscillatory (7/2 osc) and singular (7/2 sing) source terms for the equation on the acceleration ^7/2 (inertia) of y\, 
as discussed separately in Section S. 



6.2 Heuristic discussion: Time evolution 

The interplay between the two previously documented regimes of oscillatory and singular behaviors results into os- 
cillatory finite-time singularities. As a result of the nonlinearity of the restoring term (n > 1), the oscillations have 
local frequencies modulated by the amplitude of y\. We stress again that the solution yi(t; yo, to) is controlled by the 
initial condition (y , to)- In this heuristic discussion, we shall use the simplified notation y\ for y\{t\ y , to)- 

A naive and approximate way to understanding the origin of the frequency modulation is that the expression 

n-1 

7/1 defines a local frequency proportional to 1/7 1 7/1 1 2 : the local frequency of the oscillations increases 



(7 1 2/1 



In- IN 



with the amplitude of 7/1. It turns out that this naive guess is correct, as shown by the expressions ( |1 12[ ) with (|113j) of 
section [7.4.4j . We thus expect the local frequency to accelerate as the singular time t c is approached, if the amplitude 
1 7/1 (t) I grows like (t* — t)~ z (see the derivation leading to (U8), then the local period, corresponding to the distance 

between successive peaks of the oscillations, will be modulated and proportional to (t c — t)~ 



z(n-l) 
2 
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6.2.1 Case m = l 



We study 

d y i • r 11 in 

—j^T = a ~jT - 7 sign[yij|yi| , (66) 

and re-introduce the parameter a to allow us investigating the effect of its sign. 

The interplay between the l.h.s. and the first term of the r.h.s. of ( |66| ) leads to an exponentially growing trend 



^ and thus an exponentially growing typical amplitude of yi(t). If yi(t) ~ e at , both and are of the same 
order while the reversal term is of order y" ~ — e nat , showing that the oscillations will be a dominating feature of 



the solution. This is indeed what we observe in figure |6J which shows the solution of ( |66D for the parameters a = 1, 
m = l,n = 3, 7 = —10, yi(t = 0) = 1 and y 2 (t = 0) = 5. The amplitude of yi(t) grows exponentially and the 
accelerating oscillations have their frequency increasing also approximately exponentially with time, in agreement 
with our qualitative argument. 

6.2.2 Case m > 1 

Case a > and 1 < m < 2: \yi(t — ► t c )\ = \y2ft — ► t c )| — ► +°o In this regime, diverges on the approach 
of t c as an inverse power of (t c — t). The accelerating oscillations are shown in figures [7] and|| for the parameters 
m = 1.3, n = 3, a = 1, 7 = 10, y\(t = 0) = 1 and y 2 (t = 0) = 1. We observe that the envelop of yi(t) grows 
faster than exponential and approximately as (t c — where t c 4. In figure & \yi(t) \ is represented as a function 
of t c — t where t c = 4. A double logarithmic coordinate is used such that a linear envelop qualifies the power law 
divergence. The slope of the line shown on the figure gives the exponent —1.5 which is significantly different from 



the prediction —(2 — m)/{m — 1) = —2.33 given by (54) on the basis of the trend term only, i.e., by neglecting the 



reversal oscillatory term. The reversal term has the effect of "renormalizing" the exponent downward. Notice also 
that the oscillations are approximately equidistant in the variable hi(t c — t) resembling a log-periodic behavior of 
accelerating oscillations on the approach to the singularity. Here, we shall not dwell more on this regime which gives 
divergent y\ and j/2 an d concentrate rather on the rest of the paper (except for the next subsection) on the case m > 2. 

Case a < and 1 < m < 2: power law decay Equation (|2|) obeys the symmetry of scale invariance for special 
choices of the two exponents m and n. Consider indeed the following transformation where t is changed into Xt' and 



yi is changed into fiy[. Inserting these two changes of variables in (32) gives 

d 2 y[ _ rdy[ { dy\ 
dt' 2 a X m dt' l dt' 



2 d 2 y[ H m dy[ dy[ _ x , , t 

^ A = a T^,-nr\-nr\ -wvM ■ ( 67 ) 



We see that (67) is the same equation as (32) if 



for which we also have 



m — I 

n = 1 + 2- , (68) 

2 — m 



H = A ^ . (69) 



The condition (p8b holds for instance with m = 1.5 and n = 3. When the relationship (68) is true, the two equations 



(67) and (p2\) are identical and their solutions are thus also identical for the same initial conditions: y\{t) = y'\{t'). 



This implies that the solution of (32) obeys the following exact renormalization group equation in the limit of large 



times when the effect of the initial conditions have been damped out: 



Vl(t) = (Xt) , (70) 
I 1 
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where A is an arbitrary positive number and /x(A) is given by (69). Looking for a solution of the form y\ ~ t 13 , we get 



P = T - (71) 

n — 1 

This exact solution, describing the asymptotic regime t — » +oo, corresponds to the decaying regime obtained when 
a is negative and will not be further explored in the sequel which focus on the singular case a > and m > 1. 



Case m > 2: — > i c )| < +oo and \y2{t — ^ t c ) I — ► +°o In this case with a > 0, the solution of (|47|) gives a 
singularity in finite time with divergence as ~ (i c — t) -1 ^" 1-1 ). Since l/(m — 1) < 1, yi(t) remains finite with 
a singularity in finite time of the type 

m — 2 

yi(t)~y c -A(t c -t)™-i (72) 

with infinite slope but finite value Y at the critical time t c since (m — 2)/(m — 1) > 0. 

The consequence is that there can be only at most a finite number of oscillations. Indeed, since y\(t) goes to 
a finite constant, it becomes negligible compared to its first and second derivatives which both diverge close to t c . 



m-2 
m—1 



Therefore, the two first terms in ( |32[ ) dominate close to the singularity and the oscillations, which are controlled by 
the last term, finally disappear and the solution becomes a pure power law ( f72|) asymptotically close to t c . 

Figure |9] shows the solutions obtained from a numerical integration of ( |32| ) with m = 2.5 yielding the exponent 
1/3, for n = 3 and initial value y±(t = 0) = 0.02 and derivative ^-\t=o = 1/2.0 = —0.3 for two amplitudes 
7 = 10 and 7 = 1000 of the reversal term. Notice the existence of a finite number of oscillations and the upward 
divergence of the slope. As expected, the stronger the reversal term, the larger is the number of oscillations before 
the pure power law singularity sets in. The number of oscillations is very strongly controlled by the initial value of 
the slope ^-|o- F° r m = 2.5, n = 3 and initial value yi(0) = 0.02 with 7 = 10, for instance increasing the slope 
in absolute value to ^-|o = — 0.7 gives a single dip followed by a power law acceleration. Intuitively, the number of 
oscillations is controlled by the proximity of this initial starting point to the unstable fixed point (0,0), the closest to 
it, the larger is the number of oscillations. 

These properties are formalized into a systematic dynamical system approach in section 7. 

6.3 Heuristic discussion: Phase space 
6.3.1 Properties in the phase space 

Having understood the dynamics of the two elements separately (Section |5|) and with the qualitative insight provided 
by the previous examples, we pose the natural question: 

Q: Can oscillatory and/or singular dynamics persist in the presence of their interaction? 



The direct numerical integration of the equations of motion in Section |62| suggests a positive answer. In the next 
section [7[ we address this question in a formal way and construct a precise phase portrait of the overall dynamics for 
given (n, m, 7). Here we articulate the problem by identifying fundamental properties of a trajectory y(t; yo, to) in 
the phase space. First, we recall that the full dynamical equation is invariant under parity symmetry in phase space: 
(y, ^y) (-y, ~ly)- The ori g in y = (0, 0) is a fixed point: 

|y|(o,o)=0. (73) 

More precisely, it is a clockwise unstable nonlinear focus in the y phase space because the flow is divergent every- 
where: 



d d 



V • = TT-^osc = m\y 2 \ m - 1 > . (74) 
at oy2 
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Therefore, starting extremely near the origin |yo| <C 1, y(t; yo, to) undergoes a clockwise oscillation with increasing 
amplitude. 

Next, we examine the properties of the oscillations. When only the oscillatory term y 2 osc is present (i.e., rising = 
0), Section 5.1.3 showed that the amplitude and the period of the oscillations can be determined by H alone. When 
the singular term ?/2sing is added, H is no longer conserved along any trajectory but increases instead: 

j t H(y; n, m, 7 ) = (j^i) y 2s in g = M m+1 > . (75) 

The growth rate of H also increases as H increases, because a higher H corresponds to a wider range of y2 during an 
oscillation cycle, as seen from (|43|). The higher H becomes, the more effective is the impact of rising, especially once 
|2/2(*;yO)*o)| reaches 0(1). The region where |y2(t; yo, to)| < 1 corresponds to the source strips (Definition 
5.2.3| ) for the case with only the singular term ?/2sm g - 

Remark 6.3.1 

For n > 1, y(t; yo, to) nas the following characteristics: 

1. From (|46|), the frequency of oscillation increases from to oo as H increases. As a consequence, for | t/2 | <C 1 
where jlH is very small, y(t;yo,to) undergoes extremely slowly divergent oscillation with increasing fre- 
quency. 



2. From (|43|), the amplitude of the oscillations monotonically grows as H grows to oo at an increasing rate. As a 
consequence, |y(t;yo,to)| eventually diverges to oo. 



3. From (|75|), y(t; yo, to) starting from any yo approaches the origin backward in time. 

4. Therefore, |y(t;yo,to)| starting from any yo connects the origin |y(t;yo,to)| — > in backward time and 
infinity |y(t; yo, to)| oo in forward time. 

We now formulate mathematically the notion of an oscillation for the overall dynamics. When y(t; yo, to) under- 
goes a sequence of oscillation cycle in the y phase space (see for example Figure [|), y\ and 1/2 change their direction 
of motion in succession. 

Definition 6.3.2 Turn of a trajectory 

We say that y(t; yo, to) makes a turn at t = t' if the variable y\ changes its direction of motion at t = t', i.e., 

d 

^ily(t';yoA,) =^ = 0. (76) 

Each complete oscillation cycle requires two turns of y(t;yo,to)- During a time interval between two adjacent 
turns, y2 changes directions (i.e., it achieves ^1/2 = 0) if the oscillation is around the origin y = (0,0); see for 
example Figure [T[ 

Definition 6.3.3 Zero velocity curves and 

We define the two zero- velocity curves in the phase space with respect to y\ and 1/2 '■ 

F {1) = {y| ^2/1 = 0, i.e., y 2 = 0}, (77) 

F(2) ~ {y| Jt V2 = 0, Le - ™iN n " 1 = y2N m - 1 }. (78) 
F^ is nothing but the yi-axis. On the curve F^ 2 \ y\ can be expressed as a monotonic function of j/2: 

2/2 = / (2) (yi) = 7™ yiM™ -1 , (79) 
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where (yi, f^ 2 \yi)) is on F^ 2 \ An alternative way for obtaining (77) and ([78|) is to use the slope of the trajectory 
y(*;yo,*o) in phase space: 



dy_ 2 
li 



dyi 'y(*;yo'*o) 



dy2 ,dyi , 



-72/1 1 2/1 T 1 + 2/2 1 2/2 1 

y2 



(80) 



where *jg. = ±oo results in F^ and |^ = gives F^ 2 \ 



Corollary 6.3.4 Complete oscillation cycle 

Staring from a point on F^ 1 ' where y(i; y , t ) makes a turn (i.e., y € one complete oscillation cycle 

requires a set of four conditions to be satisfied in sequence: [iyi = i— > ^£/2 = i— > 4yi = i— » ^t/2 = 0]. 
Accordingly, in phase space, y(t; yo, to) cuts across [F^ 1— > F^ 2 ) F^ F^ 2 )] in succession. 



Corollary 6.3.5 Transition to non-oscillatory motion 

If y(i; yo, to) ceases to reach FW or F^ 2 \ then it can no longer oscillate around the origin. 

6.3.2 Schematic dynamics in phase space 



Having Corollaries 6.3.4 and 5.3.5 in hands, we rephrase the question in Section 5.3.1 into more specific ones: 



Ql: Under what conditions and how far does the clockwise oscillatory motion owing to t/2osc persist away from the 
origin? 

Q2: Under what conditions does the finite-time singular behavior persist, and the two singular basins resulting from 
y 2 sin g exist? 

For an intuitive grasp of issues associated with these questions, we schematically summarize in Figure [l0| the 
dynamical properties along a trajectory y(i;yo,*o) due respectively to ?/2osc and y2sing- The total velocity vector 
(dyi/dt, dyz/dt) is indicated by the arrows. The sign of the two contributing terms y2osc and fusing are given in the 
triplet {4fVi, 2/2osc : 2/2sing)- The two terms j/2osc and 2/2sing can enhance or oppose each other depending on their 
relative signs. 

Furthermore, we make the following observations. 

Remark 6.3.6 

1. The phase space is divided into six domains by F^\ F^ and y\ = 0. On F^ and F^ 2 \ the components of 
the velocity vector, Ayi and gi2/2> respectively change their sign. On y\ = and FW, Ayi and y2sing change 
sign. 

2. In the second or fourth quadrant (2/12/2 < 0) defined by F^ and y\ = 0, both t/2sing and ?/2osc have the 
same sign and hence so does 4rV2- In Figure 10, it is indicated by long thick arrow with the velocity triplet 



( 372/1) 2/2osc : 2/2sing) = ( — > ~~ : — ) in the second quadrant or (+,+ : +) in the fourth quadrant. It can be 
thought that £/2osc enhances j/2sing in a way that y(£;yo,to) flows towards one of the two (+ or — ) singular 
directions in y2 from the following reason. 

• Because y 2 sing and y 2oS c have the same sign, | j- t y 2 \ = |2/2sin g | + |2/2osc| is larger than |y 2s in g |- Therefore, if 
y(i; yo, to) remains in the second or fourth quadrant without ever leaving, it must be driven to a finite-time 
singularity with the same sign as in the case where only y 2 sing operates, but at a faster rate. Therefore, two 
singular basins inevitably exist: 
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(a) In the fourth quadrant where 1/2 > and y\ < 0: y?, — > +00 with ^7/2 > and iy\ > 0, 

(b) In the second quadrant where 2/2 < and y\ > 0: 2/2 — ^ —00 with ^2/2 < and ^yi < 0. 

'sing previously for y 2smg 



These two basins are the generalization of the two basins B^r defined previously for j/2sing only (Defini- 



tion p.2. This partially answers Q2 above. 

• The only way for y(t; y , to) to escape from the singular behavior is to leave the second or fourth quadrant 
respectively for the third or first quadrant by cutting y\ = while keeping the sign of 7/2- 

3. In the first or third quadrant {y\y2 > 0) defined by y\ = and F^ 1 ' with inside, y2osc and y2sing have 
opposite signs and the total velocity vector cannot be determined unambiguously as indicated by the velocity 
triplet (J^yi,y2osc : 2/2sing) = (+, — : +) in the first quadrant or (— ,+ : — ) in the third quadrant. When 
y(t; yo, to) enters into the first of third quadrant according to clockwise motion around the origin, it respectively 
comes from the fourth or second quadrant where ?/2smg enhances ?/2sing as indicated by plain arrowheads in 
Figure [K| Therefore, ?/2sing is dominant first. Then, the effect of y2osc gradually kicks in as |yi| increases 
towards F^ where ?/2osc balances losing as indicated by dotted line in Figure [T0|. 

• If 2/2sing remains dominant and y(t;yo,to) never reaches F^ 2 \ then jj- t yi does not change sign and yi 
keeps growing. Eventually, y2sing may completely dominate ?/2osc> and y(t; yo, to) moves quickly towards 
the terminal singularity, 

(a) In the first quadrant where y\ > and 2/2 > above F^: 2/2 - ► +00 with ^2/2 > and 4ry\ > 0, 

(b) In the third quadrant where yi < and 7/2 < below F^: 2/2 — ¥ — °o with ^7/2 < and J^ryi < 0. 
The sign of 7/2 towards the singularity is consistent for all quadrants. 

• If y(t; y, to) reaches F^ 2 \ y2osc overcomes y2sing and jj- t yi changes sign as indicated by hollow arrowhead 
in Figure [H]. Eventually y(t; y, to) has to exit the quadrant by passing FW. It follows that, if y(t; y, t ) 
reaches F^ 2 \ then it also reaches i 7 ^ 1 ) and makes a turn of an oscillation (Definition |6.3.2 ). 



In summary, the following two conditions sequentially determine whether or not y(t;yo, to) can make another turn 
starting from a turning point yo G F^ : 

1 . In the fourth or second quadrant: whether or not it reaches y\ = 0. 

2. In the first or third quadrant if it reaches y\ = 0: whether or not it reaches F^ 2 \ 

For fixed (n, m, 7), the global dynamics can be completely described by the phase portrait because this is a system 
of two-dimensional autonomous ODEs. However, this geometrical structure of the phase portrait may bifurcate as the 



value of the exponents n and m vary (Section 5.1 and 5.2). In the following section, we will examine the structure of 



the global dynamics when both elements have high nonlinearity, i.e., n > 1 and m > 2. 



7 Overall dynamics for n > 1 and m > 2 with a > 0: \yi(t — > t c )\ < +oo (except 
for isolated initial conditions) and \y2{t — > t c ) \ = +oo 



Recall from Section 5.1 for the sub-dynamical system with only the oscillatory element that the case n > 1 corre- 
sponds to highly nonlinear oscillations with a monotonically decreasing period as the amplitude of the oscillations 
increases (Figure^). From Section ^2| for the sub-dynamical system with only the singular element, the case m > 2 
corresponds to finite-time singularity with finite increment in y± and infinite increment in t/2 (Figure^). Furthermore, 
Section ^31 on the phase space of the full dynamical system showed the following results: 



1. any trajectory y(t; yo, to) starting away from the origin connects the origin in backward time and |?/2| — ► °o m 
forward time; 
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2. the oscillations may persist especially near the origin; 



3. the finite-time singular behavior should persist; 

4. theF^ 2 ) -curve is critical in determining whether or not a trajectory transits from oscillatory to singular behavior. 
As we demonstrate below, the most striking dynamical feature for the case n > 1 and m > 2 is the finite-time 



oscillatory singularity. In Section [7. 1|, we heuristically describe the global dynamics by identifying the boundaries and 
basins in the phase space using two examples. The mathematical definitions of the boundaries and basins are given 
in Section 1_2. Using the template maps, we describe the global dynamics of the boundaries in Section [7.3. 1| and 
the basins in Section [7.3.2[ Finally, we study the scale-invariant properties associated with the finite-time oscillatory 
singularity in Section [ 



7.1 Phase space description 



Two examples of phase portraits are shown in Figure 11 as a collection of trajectories with (n,m) = (3,2.5) for 



7 = 10 (Figure |1 lp-c) and 7 = 1000 (Figure |llp-f), where arrows along the individual trajectories indicate the 
direction of forward time. We observe in both examples that there are two basins (labeled by B + and B~) in the 
phase space. The superscript of individual basins correspond to the sign of the terminal direction 4fy 2 as \y 2 \ — > 00. 
The boundary between the basins is kinematically defined by the special trajectories spirally out of the origin to reach 
\yi\ — > 00 as well as 1 2/2 1 —> 00 ■ Any other trajectories result in a finite terminal value of y\ as \y 2 \ — ► 00. These 
boundary trajectories are singled out in Figures [TT| b and e (labeled by b + and b~). A typical trajectory y(i;yo,to) 
starting from yo = (—0.06, 0) are also shown in Figures [TT| c and f. 

These phase portraits confirm that oscillations indeed persist and are confined near the origin about |y| < 1 
and hence H < 1 from (35). The amplitude of the oscillations continuously grows along y(t;yo,to) m forward 
time as seen from (75). For 1 2/2 1 <C 1 starting near the origin |yo| <C 1, y(t;yo>to) is nearly vertical because of 
IJjZfel ^> l;§2/l| an d follows a constant if -curve closely (Section 5.1.3, see also Figure ^). For 1 2/2 1 no more much 



smaller than 1, H increases more efficiently by growing further in y 2 as observed in Figures |TT): and f. The period of 
the oscillations decreases continuously, because j^T{H\ n, 7) = §]j^f < along y(i; yo, to) f° r n > 1. 

Once the oscillation reaches |y| « 1 and hence H 1, the singular element y2sing works on -§3/2 more effectively. 
As a result, y(i; yo, to) starts to grow rapidly, especially when it is moving vertically in the phase space. 



Recall also that closed H contours for H > 1 are stretched out vertically for n > 1 (Section 5.1.3). Therefore, 
the oscillatory element y2osc can enhance or suppress the singular behavior of y(t; yo, to) significantly when it moves 
vertically. Such stretching effect of H is more prominent for larger 7 (compare Figures p| z to f). 

For 1 2/2 1 > 1> the dynamics is extremely singular in the second and fourth quadrants where y2osc enhances y2sing 



(Remark 6.3.6). The terminal increment of y\ is nearly zero as y 2 approaches singularity in finite time as indicated by 
the almost vertical trajectories. 

In the first and third quadrants where ^y 2 changes sign with respect to y 2oS c on F^ 2 \ the boundaries b + and b~ 
divide the phase space into B + and B~ . Any trajectory starting near b + and b~ with \y 2 \ > 1 accelerates extremely 
fast into B + or B~ , indicating that the dynamics is extremely sensitive near b + and b~ for \y 2 \ > 1. Any trajectory 
that moves away from b + or b~ for | y 2 \ > 1 does so almost vertically due to the stretched structure of If for if > 1. 
Near vertical trajectories away from b + or b~ indicate that increment of y% is finite in the first and third quadrants like 
in the second and fourth quadrants. 



7.2 Singular basins and boundary 



When the dynamics has only the trend element (Section |5.2[ ), there exist two singular basins B^ and B^ mg separated 
by a boundary 6 s ; ng determined by a collection of stagnation points where the velocity is identically zero (Definition 
5.2.1 ). In the full dynamical system, we define the two boundaries b + and 6 _ kinematically as observed in Figure 11. 
The definition of the two basins B + and B~ follow in a natural way. For simplicity and economy in notation, we use 
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± and =p to represent two distinct cases by choosing them consistently in order. For example, by "±^4 with for 
±C," we mean: i) "+^4 with — B for +C," and ii) "— A with +B for — C." Use of ± and =p is possible because the 
full dynamical equation is invariant under parity symmetry. 



Definition 7.2.1 Boundaries b + and b~ . [see Definition |5.2.l[ ] 

We define the boundary by the special trajectories that connect the origin y = (0, 0) to (+oo, +oo) and to 
(— oo, — oo) (see Figures [TTj and 12). Any trajectory y(t; yo, to) starting on yo G b^ 1 will remain on it in forward and 
backward time: 



ifyoG^, then y(t; y , t ) G b^ 

with y(-oo;y ,t ) = 



(0,0), 



for t G (— oo, oo) 

y(oo;y ,to) = (±oo,±oo) 



(81) 



Corollary 7.2.2 Asymptotic behavior ofb + and b~ . 

As |j/2 1 — * +oo, b^ asymptotically approaches where = (Figure 12). A trajectory y(t;yo,to) with 
yo G b^ 1 has the following properties. 

• It can never reach F^ 2 \ because if it does, it will have to exit the region where yiy 2 > and this leads to a 
contradiction (Remark |6.3.6 , Item 3), 

• It must stay near so that the near zero velocity 4ry2 ~ keeps the trajectory from growing rapidly to a 
singularity in a finite time. This also leads to a contradiction (Theorem [7.2. lh . 

Because the boundary is kinematically defined, we have the following theorem for the basins. 



Theorem 7.2.3 Basins B + and B~ . [see Definition |5 .2. lfl 

There exist two distinct basins B + and B~ kinematically divided by b + and b~. Any trajectory y(t; yo, to) with 
yo G B^ will remain within it and reaches a finite-time singularity, i.e., 

ifyoG^, then y (t; y , t ) € B ± for t e (-oo, t c (y ) + 1 ) 

with lim y(t;y ,to) = (0,0), 

t— > — oo 

lim y(t;yo,*o) = y c (yo), y c (yo) = (yi c (yo),±°°)- (82) 

where y c (yo) and t c (yo) are the finite-time singularity and finite-time singular interval. They depend only on the 
initial condition yo because the system is autonomous. Any other trajectory not starting in B^ initially can never 
enter into B^ by cutting across due to the uniqueness of the solution. 



As seen on the two examples in Figure [11], in the presence of both the trend and the reversal terms, the oscillatory 
behavior persists near the origin. Technically, whether or not the oscillations persist depends on the competition 
between the oscillatory and singular elements with respect to the time-scales. In Section |7.4| , we will examine this 
issue in details using scaling arguments. Here, we proceed with our discussion assuming that such oscillations do 



exist, as observed in figure 1 1 



Remark 7.2.4 

For |y| 3> 1 outside of the oscillatory region, two basins B are clearly visible (see Figures 11 and |12|): B + lies 
"above" b + , and B~ lies "below" b~ . This description is carried into the oscillatory region using the direction of the 
flow as follows. 

1. B + basin: "above" the boundaries b^, i.e., to the left of b + and to the right of 6~ with respect to the forward 
direction in the flow (see also in Figure 11). Any trajectory y(t; yo, to) with yo G B + goes to U2{t c ; yo, to) —>■ 
+oo. 
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2. B basin: "below" the boundaries b^, i.e., to the left of b and to the right of b + with respect to the forward 
direction in the flow. Any y(t; yo, to) with yo £ B~ goes to 2/2 (tc! yo> to) — ► —00. 

In other words, lies to the right of with respect to the forward direction of the flow. 
7.3 Global dynamics 

7.3.1 Dynamical properties along the boundaries 

The structure of the phase space is completely governed by the boundary 6 =t . Therefore, we first study the dynamical 
properties along b . 

Definition 7.3.1 Exit turn p ±0 as the intersection ofb^ with the y\-axis 

We define a point p ±0 as the out-most intersection of b^- with the yi-axis as shown in Figure 13. Therefore, 



y(t; p ±0 , to) makes no further turn for t > to (Definition 6.3.2 ). We call p ±0 the exit turn point. The transition from 
oscillatory to singular behaviors occurs at p ±0 . 



Definition 7.3.2 Reference trajectory y =t (t), turn points p ±k and exit time t^ k 

We define y ± (t) as the reference trajectory on b^ 1 which goes through the exit turn point p ±0 at time t = 0, i.e., 

y ± (t) =y(t;p ±0 ,0) . (83) 



In backward time, y ± (t) makes a turn by intersecting the yi-axis (Definition 5.3.2). At time t ±fc (< 0), y =t (t) makes 
the A;-th (backward) turn: 

y ± (t ±fc ) = P ±fe , (84) 

where 

p ±k = (yf k ,0) eb ± (85) 

is defined as the k-th turn points. It is located at an intersection of b and y\ axis (Figure |l3|). By construction, a 
trajectory y(t; p ±fc , t ±k ) is on the reference trajectory. Starting from p ±fc , the trajectory makes k turns before reaching 
the exit turn point p ±0 at time after a time interval — i ±fc (> 0). It does not make any more turns for t > 0. We call 
t ±h (< 0) the k-th exit time. 

Definition 7.3.3 Template map for the dynamics associated with p ±fc along b^ 

We define the template map of the dynamics along each boundary b^ 1 using the sequence of turn points p ±fc : 

■■■> P ±fc+1 > P ±fc > • • • > P ±0 • (86) 
By construction, there is no other turn points between any p ±fe + 1 and p ±fc along 6 =t . 

Remark 7.3.4 

As shown in Figure 13 , 

1. Along b^ 1 in forward time, the turn points p ±fc jump between y\ > and y\ < as y ± (t) oscillates around the 
origin: 

• along 6+: y+ (2/) < and y^ (2 ' +1) > 0, 

• along b-\ y~l (2l+l) < and y^ m > 0. 
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2. On the yi-axis, the turn points alternate between b + and b : 



yi° < Vl 1 < yi 2 < y x 3 < yf 4 < 



(87) 



Remark 7.3.5 

Three main dynamical properties associated with p ±fc are as follows: 

1. the alternating signs of yf k , i.e., (—l) k yf k < along b + and {-l) k y^ k > along b~ (Corollary |7.3.4| ); 



2. the finite number of turns Ni 



turns 



k in forward time before exiting from the oscillatory region (Corollary 



7.3.2); 



3. the exit time t ±k (Definition 



These properties are summarized in Table [j] (see also Figure 13) using the template map (>) to show the dynamics in 
forward time. 
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t^ 1 
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Table 1: Dynamical properties associated with the template map of p ±fc along the singular boundary 6* (Remarks 



7.3.5). 



We now present the dynamical properties of arbitrary trajectories along b^ 1 using the template map. To do so, we 
partition b^ into boundary segments bounded by the points p ±( - k \ 

Definition 7.3.6 Boundary segment 

We define the boundary segment A6 ± ( fe+1, * : ^ of b^ to be a segment between two adjacent turn points p ± ( fc+1 ) and 
p ±fe (both exclusive) : 



Ab ±(k+i,k) = | y ±( t ) j fo[t G ( t ±fc+i^±fc)} 



(88) 



(see also Figure |l3j). The symbol "(" and ")" in the superscript are used to indicate that both the left and right 
endpoints are outside of the interval. 



Remark 7.3.7 

1 . The semi-infinite unions of the boundary segments together with the semi-infinite unions of turn points forms 
the boundary: 

b ± = U^ =0 A6 ± ( /c+1 ' fc ) U A^-'°) Ug° =0 p ±fc , (89) 



where 



Afe ±(0,-) = { y = y ±(f) ; for t > 0} . 
2. By construction, y ± (i) makes no turn over any boundary segment A6 ± ( fc+1,fc ). 



(90) 



26 



The complete description of the dynamical properties along follows. 

Corollary 7.3.8 Dynamical properties of a trajectory y(t; yo, to) on b^ 1 

Let us consider a trajectory y(t; yo, to) on starting from a point y in the boundary segment y G Ab ±( - k+1 ' k ^ 
at an arbitrary time to- It reaches the turn point p ±fc in forward time after a time interval — At^ k+1 ' k \> 0) without 
making any turn, where 

f t(fc+l)_ f tfc <Ai ±(fc+i,fc) <0) (91) 

because t ±k - t ± ( fc+1 ) is the total time of flight over Ab ±( - k+1 ' k \ Accordingly, y(t; yo, to) with y G Ab ±( - k+1 ' k ^ has 
the same dynamical properties as y^ 1 (t) starting from p ±fe (Table [I]), except that the exit time is extended from t ±k to 
t ±K + Atj v ' ; . It can also be expressed using the reference trajectory 

y(t; y , to) = y ± (t-t + t ±k + Atf {k+1 > k) ) , (92) 

because y(t ; y , to) = y±(t ±fe + At± (fc+1 ' fe) ). 



7.3.2 Dynamical properties in the basins 

As the dynamical properties along the boundaries b^ can be described by a template map of turn points p ±k , the 
dynamical properties in the basin can also described by a template map of turn segments as follows. 

Definition 7.3.9 Turn segments Ae ( - ±( - k+1 ^ k ^ 

Any trajectory y(t;yo,^o) can make a turn only on the j/i-axis (Definition 6.3.2 ). We define a turn segment 



Ae (±(fc+i),Tfc) as being bounded by two adjacent turn points p ± ( fc+1 ) and p Tk on the j/i-axis (Figure 13): 
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The notational convention for the superscript of the turn segment Ae^ k+1 ''^ k ' is that the left and right indices, ±(fc+ 
1) and =F&, respectively correspond to the superscript of the turn points, p ± ( fc+1 ) and p Tk , which are respectively at 



the left and right ends of the segment with respect to forward direction of the flow (see Remark 7.2.4). The symbols 
"(" and ")" in the superscript mean that these intersection points are not included in the segment. 



Remark 7.3.10 

In comparison with Remark 7.3.4 : 

1. Along the flow in forward time, the turn segments jump between 2/1 < and 2/1 > as a trajectory in B^ 
oscillates around the origin: 

• in B+: 2/1 < for Ae(+( 2,+1 )'-( 2 9> and 2/1 > for Ae^ 2 ').-^ 1 ) ; 

• mB~: 2/1 < Ofor Ae^^)^ 2 '- 1 )) and 2/1 > for Ae (-( 2m )<+( 2/ ); 
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Moreover, the yi-axis consists of the union of all intersection points and segments; 

yi -axis = Ae^ - 00 ) U p+° U Ae^'+V U p 1 U Ae^ 2 '" 1 ) U p+ 2 U 

. . . U (0, 0) U . . . U p 2 U Ae ( ~ 2 ' +1) U Ae (+1 '- Q) U p~° U Ae ( -°' +oo) . (93) 

Note that the superscript of Ae( ± ( fe+1 ) ,=Ffc ) for y\ < is not in sequence with the subscript of p ± ( fc+1 ) as in the 
case for y\ > 0. This is because p^ -1-1 ) and p^ k are located at the right and left ends of Ae < - ±( -' c+1 - ) ' =Ffc ^ with 



respect to the forward direction of the flow (Definition 7.3.9), but geometrically at the left and right ends of the 



segment on yi-axis (Figure 13). 



To describe the dynamics in the basins, we first define the oscillatory source near the origin and separate it from the 
singular region outside. We show that the transition from the oscillatory to the singular behavior in B occurs at 
the exit turn segment, associated with the fact that the transition along occurs at the exit turn point (Definitions 
[7.3.1| and 7.3.2). The global dynamics in B^ is structured into two regimes separated by the exit turn segment which 



determines the singular behavior in forward time and the oscillatory behavior in backward time, as follows. 

Definition 7.3.11 Exit turn segments Ae^ 1 ^ ' G B ± 
We call Ae( ±1,T °) G B^ the exit turn segments. 

Definition 7.3.12 Oscillatory Source S 

The exit turn segments and the out-most boundary segments, along with the turn points at the end of these 
segments, formally define the oscillatory source region S: 

S = {y| region surrounded by p _1 , Ab^^'°\ p _0 , Ae( +1,_0 ), 

p+\ A6+( 1 -°), p+°, Ae^ )} (94) 



(see Figure 13; as well as Definition 5.2.3 for the case with only the singular element). 



Corollary 7.3.13 Transition from oscillatory to singular behavior 

In the sequel, we note p( ±1 ' = F0) f or short to include one of the four points p T °, p^ 1 . Let us consider a trajectory 
y(t; yo, 0) starting from a point yo = p( ±1,T °) on the exit turn segment Ae^ 1 '^ ^ G B^. In forward time, y(t; yo, 0) 
will make only one turn at a point in Ae^ ±0 ' T °°^ but never completes an oscillation cycle (Corollary |6.3.4 ). Therefore, 
the exit segment Ae^ ±1,=F0 ^ defines the transition from oscillatory to singular behavior. 

Corollary 7.3.14 Dynamics outside S in the singular region. 

Let us consider a trajectory y(t; yo, 0) starting from the exit turn segment in forward time with yo = p^ 1 ^ ) g 
Ae( ±1,=F °) in B^. After making the final turn in Ae^ ±0 ' TOO \ it reaches the corresponding singularity: 



y c (p 



(±1,=f0)n 



y^p^ ));?^ 1 '^) £B ± (95) 



which depends only on the initial condition p( ±1 '^°) ( se e Theorem |7.2.3| ). Here t c (p( =tl ' =F °)) is the finite singular 
time. 

By Definition 7.3.9j , the left and right end points of Ae ( - ±1,T °' 1 are next to p^ 1 and p T °. Therefore, the terminal 



value of y\ at the singularity ranges over: 

yie(p (±1 ' T0) ) = <±oo, T oo> , (96) 

where (a, b) denotes that it is a or b asymptotically if p( ±1 ' = F0) j s respectively at the left or right end point of Ae^ 1 ^ ' 
with respect to the forward direction of the flow. 
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Remark 7.3.15 

Two main dynamical properties associated with a point p( ±1 '^o) on the exit turn segment in the singular region outside 
S are as follows: 

1. finite singular time, t c (p( ±1,=F0 )); 

2. terminal y\ value, 2/i c (p^ ±1,T0 ^)- 



Corollary 7.3.16 Dynamics inside S in the oscillatory region. 

Let us consider a trajectory y(t; yo, to) starting from the turn segment Ae^ ±k ' T ^ k ^ 1 ^ G B^ in forward time . In- 
cluding the starting point as the first turn, the trajectory makes the Z-th turn (1 < I < k) at a point in /\e( +k ~ l+1 ~( k ~ 1 ^ 
and the sign of y\ alternates between + and — at each turn (Remark 7.3.1C ). The trajectory reaches p( =tl ' = F°) f an 



exit turn segment Ae^ ±1 ' T °' 1 to make the k-th and final turn at time to + At^ k ' T ^ k ^\ with 



(97) 



where t ±k are defined by (Q) (see also table [j] (a, 6)t denotes that it is a or 6 asymptotically if yo is respectively at 
the left or right end of AeW^ 1 ''^' with respect to the forward direction of the flow, like (a, b) in ( ^6| ) for y\ c . 



Definition 7.3.17 Template map for the dynamics associated with Ae^ k+l ^^ k 

We define the template map of the dynamics in using the sequence of intersection segments Ae^*" 1 " 1 ^*) 
on the yi-axis: 

... >Ae (±(*+l), T fc) > Ae (±fc,T(fc-l) >>>>>Ae (±l,TO) >Ae (±0,Too) (9g) 

By construction, there is no other turn segments between Ae^ ± ^ k+1 ^' Tk and Ae^ ±k ' T ^ k ^ 1 ^ in B^ (compare with 
Definition [7.3.31). 



Remark 7.3.18 

Three main dynamical properties associated with a point yo G Ae^ =t ^ +1 ^' =Ffc ^ in the oscillatory source S are summa- 
rized in table |2| (compare with Remark 7.3.5[ ). They comprise: 



1. the signs of y\ (Definition [7. 3.9| ); 

2. the finite number of turns Ntum to transit from S into the final region where the monotonous singular behavior 
occurs (Corollary [7.3. 16| ); 

3. the exit time interval Ati ±(fc+1) ' Tfe) tQ reach Ae (±i, T o) (Corollary gjOg ). 

We now present the dynamical properties of arbitrary trajectories in B^ using the template map. To do so, we partition 
B^ into sub-basins limited by the segments Ae( ± ^ +1 ^' =Ffc ^. 

Definition 7.3.19 Sub-basins AB ± ^ k+1 '^ 

We define a sub-basin to be a piece of the basin B^ limited by two adjacent turn segments Ae^ k+1 '^ k > and 
Ae (±fe,=F(fe-i)) as follows (see Figure ||). 

A ^±(fc+i,fc) _ |y| re g lon exclusively surrounded by (99) 

Ae (±(fc+l),Tk) ±fc A6 T(fc,(*-l)) p ±(*+l) 
Ae (±fc, T (fc-l)) 5 pTfc^ A6 ±((ft+1),*) p =F(fc-l)\ . 
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Vic 


yic(p^W) 


S 


tc 


t c (p(±^u)) 



Table 2: Dynamical properties associated with the template map of Ae^ k+1 ^' Tk ^ in the oscillatory (in S) and singular 



(outside S) domains. These properties characterize the dynamics in forward time (see also Figure |13| and Table |lj) 
using the template map (>). 



Remark 7.3.20 

In comparison to Remark 7.3.7 : 



1. The semi-infinite unions of the sub-basins together with the semi-infinite unions of turn segments reconstruct 



the complete basin (see Figure 13): 

B ± = U^ =0 AB^ k+1 ^ U AS ± (°'°°) u£L Ae< ±fe+1 ' fc) , (100) 

where 

region exclusively suiTounded by (101) 
A&+(°<-\ p +0 , Ab+^\ p+\ Ae( +1 '-°), p-°, Afe-(°'-)} 



where Afe^ '") are defined by 
2. By construction, no trajectory makes any turn in AB^ k+1 ' k \ 

The complete description of the dynamical properties in can now be given completely. 

Corollary 7.3.21 Dynamical properties of a trajectory y(t; yo, to) m B 

Let us consider a trajectory y(t; yo, to) in B^ starting from a point yo inside the sub-basin yo G AB^ k+1 ' k ^ at an 
arbitrary time to- It reaches the turn segment Ae^*'^^ -1 ^ in forward time after a time interval — At^ k+1 ' k \> 0) 
without making any turn, where 

t ±(k + i) _ t ±k < At ±(*+l,*) < o f (102 ) 

because t ±k - t ±( - k+1 \> 0) is the time of flight of a trajectory along A& ± (( fc+1 )> fc ) which borders AB ±( - k+1 ' k \ 
Accordingly, y(t;yo,to) with yo € AB^ k+1,k ' has the same dynamical properties as Ae^*'^* -1 " (Table ||), 
except that the exit time is extended from Ati ±k ' T ^ k 1 ^ to At e ±k ' T ^ k ^ + Atg^ k+1 (compare with Corollary 



7.3.8). 



7.4 Scaling laws 

7.4.1 Dynamical properties on the yi-axis 

Because the system is autonomous, each trajectory y(t; yo, to) is determined uniquely by its initial condition yo- We 
have shown above that the template maps defined on the yi-axis completely describe the dynamical properties of the 
dynamical system. It is thus convenient to summarize these dynamical properties as a function of the initial condition 
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Yo = (yiO) O) taken on the yi-axis. There properties are quantified by the finite singular time t c , the total number of 
turn N tU rn, and the finite terminal value y\ c (see Remarks 7.3.5| and 7.3. 18| ). 



Figure [14] shows these dynamical properties as a function of y\Q for (n,m) = (3,2.5) with 7 = 10 and 1000, 
determined by direct numerical integration of the equations of motion, using a fifth-order Runge-Kutta integration 
scheme with adjustable time step. Each discontinuity of N turn and y lc as a function of t/ 10 occurs at a turn point p ±fc 
which separates two turn segments, one in B + and the other in B~ (Figure Hand Remark [73l0| ). 

Ntum is directly associated with the oscillatory behavior of the dynamics. Note that it exhibits a staircase structure, 
being constant for any point in Ae^^ k+1 ^' Tk ^ (Table ^). 

In contrast, y\ c is a property more directly associated with the singular behavior and takes identical values for 
any point yo € /\ e ±k,T{k-i) which is mapped to yi c (p^ =tl ' =F °') = (±00, =1=00) after k turns (Corollary [7.3.14 ) with 



yic(p ±fc ) = ±00. Moreover, given a y\ c , each turn segment Ae( =t ( fc+1 )' =F ' C ) has a unique point y ± ( fc+1 )' =Ffc ) (Corollary 



The critical or singular time t c is function both of the oscillatory and the singular terms: 

i c (yo)=t c (p (±1 ' TO) )-A4 ± ( fc+1 )^), (103) 

where yo G /\ e ( ±( y k + l )^ k ) . However, the explosive singular time scale i c (p( ±1,=F °)) is generally much shorter than 
the slow oscillation time scale L 
is dominated by A^ ±(fc+1) ' Tfc) . 



the slow oscillation time scale Aii ' (see figure 0)). Hence, the total time t c needed to reach the singularity 



7.4.2 Definition and mechanism 



Our analysis up-to-now has demonstrated a hierarchical organization of the spiraling trajectories diverging away from 
the origin in phase space, as shown in figure [13|. Figure 14 quantifies this hierarchical organization by showing the 
dependence of the critical time t c , the number N tuIIl of rotations of the spiraling trajectory in phase space and the 
final value y\ c , as a function of the initial value of yio = yi(io) on the yi-axis. The two former quantities diverge as 
power laws of y\Q with negative exponents as y\Q — > 0. The last quantity exhibits a "local fractal" structure around 
the origin which reflects the nested spiral structure of the two basins B + and B~ around the origin S, and the fact that 
each turn segment Ae^^ 1 )^*) shares the same singular dynamical properties as Ae( ±1,T °\ Accordingly, iVtum is 
k + 1 for y 10 € Ae^ k+1 ^ 



Figure 15 make these statements more quantitative by showing the log-log plots of i c (yio)> of At c (yio) (defined 
as the increment of t c over one turn of the spiral starting from a given initial point), of iV turn (yio) and of the increment 
Ayi(yio) over one turn of the spiral, as a function of (yio)- The observed straight lines qualify power laws. In order 
to get accurate and reliable estimations of these dependences and of the exponents defined below, we have integrated 
the dynamical equations using a fifth-order Runge-Kutta integration scheme with adjustable time step. 

Figure [1^ shows that the exponents (slopes of the log-log plots) are essentially identical for 7 = 10, 100 and 1000, 
indicating that the scaling properties depend only on the exponents (n, m) and are independent of 7. 



These different power laws correspond to the linear behaviors shown in figure |15| and can be represented as 
follows. 



-Wturn 
tc 

At c 



V10 
V10 



where a > , 
where b > , 
where c > , 
where d > , 



(104) 
(105) 
(106) 
(107) 



t c (p+ fe )|,t c = i c (p +fe )>andiV I 



where j/iq = y\{ta). Specifically, the definitions are (see figure 16): Ay% = \yiQ k+2 ^ — yfo\, At c = |£ c (p + ( fc+2 ^) 



tarn 
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The self-similar behavior and power laws occur because there is a countable infinite number of /\ e ^( k + 1 < k ) within 
an extremely slowly-divergent oscillatory source region S which reach the pre-exit basin segment Ae( ±1,=F °) after k 
mapping (therefore after about t ± ( fc+1 ) exit time interval with N tnrn = k + 1 forward turns). 

Note that the self-similar behavior close to the origin is governed by the singular boundaries b^. The choice of 
the z/i-axis to define the segments has been made for the simplicity associated with counting number of forward turns. 
In theory, ?/i-axis can be replaced with any judiciously chosen curve/lines. Here, we choose the yi-axis to define the 
segments, so that the self-similar properties are described as a function of the initial position y w while its velocity y 2 
is set to zero. 



7.4.3 Scaling relations from self-consistency 



Eliminating y\Q between (105) and (107) gives 



(108) 



Since At c is nothing but the difference At 
on the function t c (N tmn ) 

Upturn + AiV turn 



t c (Ntum + 1) — ^c(-^Vturn). ( |108 ) gives the discrete difference equation 

tc(Ntum) J 



AN t} 



(109) 



where AN tUTn = N tmn + 1-Ni 



turn 



1. The left-hand-side of (11091) provides a discrete difference representation of 

i d 



the derivative dt c /dN tUTn . Equation ( |109[ ) can then be integrated formally to give N turn ~ t c b , which is valid for 
d < b. Comparing with the relation between /A/turn and t c obtained by eliminating y\o between ( |104{ ) and d!05| ), i.e., 

a * ' * * 

Ntum ~ tc , we get the scaling relation 

a = b-d. (110) 

Since a > 0, the condition d < b is automatically satisfied. 

Similarly, Ay\ = yio(N tmn + 1) — yio(N tuin ) oc y\ Q according to definition ( |106| ). This gives the differential 
equation dyio / dN tmn ~ y± , whose solution is A^um ~ vau d for c > 1. Comparing with the definition 



( |1()4| ), we get the second scaling relation 

a = c 

Since a > 0, the condition c > 1 is automatically verified. 



1 . 



(Ill) 



The scaling relations (111) and ( 1 1 1 ) are the only two that can be extracted. This shows that the exponents defined 



in ( 104 ■ 107) are not independent: out of the four exponents a, b, c, d, only two of them are independent. 



7.4.4 Determination of the critical exponents 



To go further, we use the insight provided by sections a to 7.3 



Deep in the spiral structure described in the previous section and depicted in figure 13, one full rotation would 



close on itself in the absence of the trend term and would conserve exactly the Hamiltonian H defined in (35). In 



this case, we know that one full rotation takes a time equal to the period T{H) given by (44). In the presence of the 
trend term, one full rotation is not closed but the failure to close is very small especially so very close to the origin. 
Actually, the failure of the trajectory to close is quantified by the variable Ay\ introduced in the previous section and 
used in ( JTO^ ). 

The approximation we are going to use is that the value of the period T(H) without the trend term gives the 
time At c needed to make one full rotation (notwithstanding the fact that it does not exactly close on itself). This 
is essentially an adiabatic approximation in which the Hamiltonian H and the period T(H) are assumed to vary 
sufficiently slowly so that the local period of rotation follows adiabatically the variation of the Hamiltonian H. 
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Equation (14) gives T(H) oc i? 2 (™+ 1 ) . From the normalization (36), the amplitude maximum of y\ is proportional 



to H "+ 1 . Putting these two relationship together gives 

At, - 



which, by comparison with ( |107| ) gives 



d 



1 2/io I 



n 



(112) 



(113) 



We need a second equation to determine completely the other exponents a, b and c. It is provided by dT/dt, 
expressed as dT/dt = (dT/dH) x (dH/dt), where dT/dH is obtained from @ and dH/dt is given by ©. 
Estimating dT/dH from (|^) is consistent with the above approximation in which a full rotation along the spiral 
takes the same time as the corresponding closed orbit in absence of the trend term. Expressing dH/dt using (75) 



involves another approximation, which is similar in spirit to a mean-field approximation corresponding to average 
out the effect of the trend term over one full rotation. In so doing, we average out the subtle positive and negative 
interferences between the reversion and trend terms depicted in figure [Kj Furthermore, consistent again with the 
above approximation in which a full rotation along the spiral takes the same time as the corresponding closed orbit in 
absence of the trend term, we replace dT/dt by dAt c /dt c and obtain 



dAt t 

dt c . 



3n+l 



H 



12/2 



m+1 



Ate 1 



3n + l 

n-1 



(ra+l)(m + l) 

ml 2 



(114) 



where the dependence in At c in the last expression of the right-hand-side is derived by replacing H by its dependence 
as a function of T (by inverting T(H) given by (0)) and by identifying T and At c . In the last expression, we have 

r— , n+1 

replaced the dependence in 1/2 by the dependence in y\ by using the normalization (pq), leading to 1 2/2 1 ~ \y\\ 2 . 
Taking the derivative of (|108|) with respect to t c provides another estimation of 737^, and replacing At c in (114) by 



its dependence as a function of y w as defined by (107) gives finally: 



a = — (n + l)(m + 1) 



-(3n + l) 



(115) 



Figure 16 presents a comparison of the theoretical predictions ( 1 1C ), (1 1 1 ), ( 1 13 ), (UTS) for the exponents a, 6, c, d 



defined by (|104|)-(|107j) with an estimation obtained from the direct numerical simulation of the dynamical equations. 
The lines are the theoretical predictions of the exponents a, b, c, d for m = 2.5 (solid line), m = 2.75 (dotted line) and 
m = 3 (dashed line) as function of the exponent n. The symbols correspond to the exponents obtained by numerical 
simulation for 7 = 10 (square), 7 = 100 (diamond) and 7 = 1000 (crosses). The agreement is validated to within 
numerical accuracy. We verify also the independence of the exponents a, b, c, d with respect to the amplitude 7 of the 
reversal term. 



7.4.5 Time-dependent expression of the envelop of y\ (t; y , to) 



We have seen in section p.2.2| that, after exiting from the spiral in phase space, the dynamics becomes completely 
controlled by the trend term, while the reversal term responsible for the oscillations becomes negligible. This leads to 



the asymptotic solution close to t c given by expression (|72|), which we rewrite here for the sake of comparison: 



yi(t) ^ yie - A(t c - ty 



outside the oscillatory regime 



(116) 



where A is an amplitude. yi(t) has an infinite slope but a finite value y\ c at the critical time t c since < (m — 
2)/(m — 1) < 1. The dependence of this finite critical value y\ (t c ) = y\ c as a function of y\o is shown in figure [l4|. 
In the oscillatory regime, we can also obtain the growth of the amplitude of y±(t) by combining some of the 

Since Ayi 



previous scaling laws ( JlQ4j - |l 07[ )- Indeed, taking the ratio of ( |106| ) and ( |107| ) yields Ayi/At, 



2/io 
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corresponds to the growth of the local amplitude A yi of the oscillations of yi(t) due to the trend term over one turn 
of the spiral in phase space, this turn lasting At c , we identify this scaling law with the equation for the growth rate of 
the amplitude A yi of y\ in this oscillatory regime: 



dA 



yi 



dl 



J^c+d 



Its solution is of the form 



A yi (t) 



B 



(t* - ty/ b 



within the oscillatory regime 



(117) 



(118) 



where B is another amplitude. We have used the scaling relations (11C) and (111) leading to c + d — 1 = b. The 



time t* is a constant of integration such that B/(t*)t = A yi (to), which can be interpreted as an apparent or "ghost" 



critical time, t* has no reason to be equal to t c , in particular since the extrapolation of ( |118| ) too close to t* would 
predict a divergence of y\(t). The dynamical origin of the difference between t* and t c comes from the fact that t* is 



determined by the oscillatory regime while t c is the sum (103) of two contributions, one from the oscillatory regime 
and the other from the singular regime. 



The prediction (118) is verified accurately from our direct numerical integration of the equations of motion, as 
shown in figure 17. To get this figure, we rewrite ( |118| ) as 



t = t* 



(3[A yi (t)} 



(119) 



where [3 = B~ b . Note that ( |119D has the same power law dependence on y\ as (105). The reason lies in the fact that, 
for yio very close to the origin, the duration of the oscillatory regime is overwhelming that of the singular regime. As 

should and are the same. 



a consequence, the two power laws ( |1 19] ) and ( |105| ) defined in the asymptotic limit yio 

We compare this expression with the numerical simulation using \yf k \(t +k ) for A yi (t). The triangles are the 
data) (t +k , \yi k \) which are fitted to (119) to get t* and (3. The exponent b is fixed to its theoretical value given by 



( 1 10|jl 13|J1 15 ). As can be seen from the two panels, t +k as a function of \yf | is a straight line as predicted by (|TT 



with a very good precision. The first panel shows the whole calculated range. The second panel shows a magnification 



close to the exit point of the oscillatory regime. The different straight lines correspond to fits of the data with (119) 
over different intervals. We obtain respectively 



t* 


= 0.3857 and (3 = 


0.3381 for the interval k 


e 


[i- 


301]; 


t* 


= 0.6665 and (3 = 


0.3381 for the interval k 


e 


[61 - 


-» 301]; 


t* 


= 1.9577 and (3 = 


0.3380 for the interval k 


e 


[241 


-» 301]; 


t* 


= 2.4330 and (3 = 


0.3380 for the interval k 


e 


[271 


-> 301]; 


t* 


= 3.5695 and (3 = 


0.3379 for the interval k 


e 


[286 


-> 301]. 



These five fits performed increasingly deeper within the oscillatory regime exhibit a good stability for the determi- 
nation of the slope parameter (3 = 0.3380 ± 0.0001 but an alarmingly strong variation of the "ghost" critical time 
t*. Essentially, we must conclude that it is impossible to determine t* with any reasonable accuracy. There are two 
reasons for this. First, as the next section 7.4.6| shows, there is a very slow shift or cross-over within the oscillatory 
regime to the final power law singularity at t c . This cross-over corresponds to yi(t) starting off from its initial value 
with an amplitude growing according to (118) and then crossing over to ( |116| ) close to t c . Second, while we hoped 
that the asymptotic behavior ( |118| , |119 ) would become stable deeper within the oscillatory regime, this is counteracted 



by larger distances of \yt k \(t +k ) from the origin, which make the determination of t* more unstable. These diffi- 
culties are similar to but stronger than the well-known problems of the accurate determination of critical transition 
parameters. 
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Combining (107) with (118) gives the time dependence of the local period At c of the oscillation in the oscillatory 
regime: 

At c ~(t*-t)i. (120) 



Notice that the local period of the oscillation is not constant but shrinks progressively to zero as time increases. This 



is qualitatively reminiscent to the log-periodic oscillations discussed in the introduction Q49p . Quantitatively, there is a 
difference which can be characterized as follows. A log-periodic oscillation associated with discrete scale invariance 
(DSI) corresponds to characteristic time scales t ±k (for instance, the times of the local maxima of the oscillations) 
such that 



t ±(k+i) _ t ±k „ 



1/X k 



where k > is an integer and A > 1 is a prefered scaling ratio of DSI. Expression (121 ) predicts that t* — t 



i.e., goes to zero exponentially as a function of the index k. In contrast, expression (|120|) can be rewritten 



t ±(k+i) 



+ ±k 
t ~ 



(t*-t 



±k 



(121) 
(122) 



Expression (122) predicts that t* — t ±k ~ {k* — A;) 1 /*- 1 "), i.e., the local period goes to zero in a finite number of 
turns. The log-periodic result is obtained as the limit d/b — ► 1~, which is reached for n — » oo and m — ► 2. We 
thus have a dynamical theory that provides a mechanism for quasi-log-periodic oscillations with, in addition, a finite 
number of them due to the cross-over to the non-oscillatory regime. A similar almost log-periodic but nevertheless 
distinctly different frequency modulation has recently been observed numerically on the average of the logistic map 
variable close to a tangent bifurcation associated with deterministic intermittency of type I These "power law 
periodicity" are originated in the mechanism of reinjection of the iterates on the channel of near-periodic behavior. 



7.4.6 Deviation from scaling 



Figure |15| and even more so figure |16| have been constructed in the "scaling" regime in which the two approximations 
used above hold strongly. 

1. Adiabatic approximation: the value of the period T(H) calculated without the trend term gives the time At c 
needed to make one full rotation (notwithstanding the fact that it does not exactly close on itself). Equivalently, 
the Hamiltonian H and the period T(H) are assumed to vary sufficiently slowly so that the local period of 
rotation follows adiabatically the variation of the Hamiltonian H. 

2. Mean-field approximation: we average out the effect of the trend term over one full rotation. In so doing, 
we average out the subtle positive and negative interferences between the reversion and trend terms depicted in 



figure 10. This allowed us to estimate dH/dt using (75). 



The theoretical predictions ( p0| ), ( |TT1| ), flTDl ), < JTT^ ) for the exponents a, b, c, d defined by (JTCM|)-(|T07|) obtained 
using these two approximations have been found to be in very good agreement with direct numerical estimations, as 



shown in figure 16 



However, this agreement is obtained only within a "scaling" regime, sufficiently close to the origin in phase space, 
i.e., such that the Hamiltonian H of the oscillations grows sufficiently slowly. To quantify the concept of a "scaling" 
regime, figure [T^ represents the terminal critical value yi(t c ) = y\ c as a function of initial value yi(to) = 2/io m 
intervals such that perfect self-similarity can be checked readily. The top panel presents a synopsis by showing the 
oscillations of y\ c as a function of y\o for the first 300 turn segments Ae( ± ( fc+1 )' =Ffc ) (see definition 7.3.9). If self- 



similarity was true, the three following panels in figure 18 should be essentially identical because they show exactly 
the same number of oscillations. However, it is clear that there is a systematic drift, which is fast at first (second panel 
for the first 50 turn segments Ae( ± ( A:+1 )' =Ffc ) and then slows down deep within the spiral structure for the third and 



fourth panel (see figure 13 for definitions). 
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Figure 19 exhibits almost perfect self-similarity by plotting the terminal critical value yi(t c ) = y\ c as a function 
of initial value yi(fo) = 2/io from the 240th to the 260th turn segments Ae^ ±( - k+1 ^' Tk ^ (top panel), from the 260th to 
the 280th turn segments Ae^ k+1 ^ Tk ^ (middle panel) and from the 280th to the 300th turn segments Ae( ±( - k+ V>^ 
(bottom panel). The two vertical lines provide a guide to the eye to verify the almost perfect self-similarity This is the 
regime and beyond where the scaling laws reported in figures 15 and [if] hold. The cross-over to this scaling regime is 
presented visually in figure^ which compares exactly twenty oscillations of the terminal critical value yi(t c ) = y± c 
as a function of initial value yi(to) = yio in different intervals, from close to the exit point (first top panel) to deep 
within the spiral structure (bottom panel). The two vertical lines point out the deviations from the close-to-asymptotic 
regime of the bottom panel reached from the 280th to the 300th turn segments Ae^ ± ^ fe+1 ^ ,=Ffc ^ as the dynamics gets 
closer and closer to the exit point (from the third to the first panel). 

It is probably possible to improve upon the scaling theory offered in section 7.4.4 and go beyond the adiabatic 
mean-field approximation in order to describe the cross-over regime between the asymptotic scaling regime and the 
exit of the spiral structure. This should be done by quantifying the relative reinforcing and opposing effects of the 
trend and reversal terms within each turn, as shown in figure [HI 



8 Concluding remarks 

We have introduced a second-order ordinary differential equation describing an oscillator exploding in finite time 
whose dynamics results from the interplay between a nonlinear negative viscosity and a nonlinear reversal term. This 
system provides a simplified description of stock market prices, population dynamics and material rupture, in regimes 
where the growth rates are an increasing function of the price, population or stress, respectively, in the presence of 
important negative feedback effects. Our approach using dynamical system theory has shown that the trajectories can 
be understood in details from the spiraling structure of two sets of specials curves in phase space linking the origin to 
points at infinity. 

The message of our work is threefold. First, it is important to take into account the delayed response of dynamical 
variables to past states, leading technically to the presence of a second (or higher) order differential equation. This 
inertia is essential for the generation of oscillations resulting from overshooting. Second, positive nonlinear feedback 
is a general and ubiquitous mechanism for generating accelerated super-exponential growth ending in finite-time 
singularities. Third, reversal, recovery or healing mechanisms with nonlinear threshold-like behavior, together with 
inertia, ensure the existence of overshooting and thus of oscillations controlled by the amplitude of the variable. Our 
two-dimensional nonlinear dynamical system opens the road to a re-examination of many systems which contain these 
elements but which have been linearized, thus missing the novel phenomenology unraveled here. Our study suggests 
that super-exponential growth modulated by amplitude-dependent oscillations may be a general feature of complex 
systems, such as financial markets, population and heterogeneous materials. We believe that this work provides a 
novel framework to model these systems and to discover new precursory indicators or patterns of "rupture". This 
seems to offer a generalization to the reported log-periodic critical oscillations previously reported for these systems 
(see [|^] and references therein). 

The study presented here can be enriched in many ways. For several applications, three obvious missing ingredi- 
ents are additive noise, stochastic reversal to the fixed point (yi = 0, = 0) and saturation effects with reinjection 
limiting the growth before reaching the singularity. The addition of these terms will naturally lead to intermittent 
oscillatory structures (deterministic, stochastic, or both resulting from the interplay between deterministic chaos and 



noise), similar to those documented as log-periodic power laws in financial bubbles \ pl[ y.l\ [28| ] as well as in the 
population dynamics [BTJ and in rupture [0, E: 
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Figure 1 : Normalized model for the oscillatory term: phase space (left panels) and time series (right panels) for a,b) 
n = 0.5; c,d) n = 1; e,f) n = 3; and g,h) n = 15. The continuous (resp. dashed) line is yi (resp. y{). 
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Figure 2: Phase portrait (left panels) and period T of the oscillations (right panel) for the oscillatory term with 7 = 1 
for all, so that H = 1 corresponds to the normalized model (Figure |l|): a,b) n = 0.5; c,d) n = 1; e,f) n = 3; and g,h) 
n = 15. The period of the oscillation is on the abscissa as a function of the maximum equal to {2H)i of yi on the 
ordinate. 
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Figure 3: Normalized model for the singular term: phase space trajectories (left panels), critical time t c (center 
panels) as a function of the initial value 3/2,0 an d time series (right panels) for a-c) m = 1.5; d-f) m = 2; g-h) 
m = 2.5. Thick and thin lines correspond to the pair of normalized orbits for t/2 > and ?/2 < 0, respectively, with 
initial condition (1/1,0, 2/2,0) = (0, ±0.6). The continuous (resp. dashed) line is 2/2 (resp. y\). 
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Figure 5: Two singular basins and the boundary between them for the singular term with (m, a) = (2.5, 1) (see also 
Figures ||g and fj; for the normalized model) in the y phase space. The reference orbits are labeled by 0. 
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Figure 7: Solution of (32) for the parameters m = 1.3, n = 3, a = 1, 7 = 10 and 1/20 = dy\/dt\ t= o = 1. The 
envelop of yi(t) grows faster than exponential and approximately as (t c — t) -1,5 where t c 4. 
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a=l 7^10 m=1.3 n=3 




Figure 8: Same data as in figure [7]: the absolute value |yi(i)| is shown as a function of t c — t where t c = 4 in log-log 
coordinates, such that a linear envelop qualifies the power law divergence (t c — i) -1,5 . The slope of the line is —1.5. 
Notice also that the oscillations are approximately equidistant in the variable ln(i c — t) resembling a log-periodic 
behavior of accelerating oscillations on the approach to the singularity. 
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Figure 9: Solutions obtained from a numerical integration of ( p2| ) with m = 2.5 yielding the exponent ^| = 1/3 for 

m — 2 

the terminal singular behavior of y\ ~ y\ c — A(t c — t) m ~ 1 close to t c , for n = 3, a = 1 and initial value yio = 0.02 
and derivative 1/20 = ^fr|o = —0.3 and for two amplitudes 7 = 10 and 7 = 1000 of the reversal term. 
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Figure 10: Schematic velocity field indicated by arrows. The phase space is divided into six regions by FW, 
and yi = 0. The effect of the reversions and trend terms in each region are shown by the sign of y 2oS c and losing 
in (^yi,y 2 osc : y 2s in g )- 2/2osc and y 2sing may enhance each other (long thick arrows); y 2sing dominates y 2osc (plain 
arrowhead); y 2oS c dominates y 2s in g (hollow arrowhead). On F^ 2 \ they balance (dotted line). 
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Figure 11: Global dynamics in the phase space with (n, m) = (3, 2.5) for 7 = 10 (a-c) and 7 = 1000 (d-f): a,d) 
phase portrait as a collection of trajectories; b,e) singular boundary 6 =t , and c,f) a trajectory starting yo = (—0.06, 0) 
with contours of H (dotted lines) and G (dashed lines) for reference (see also Figures ^ and ||, respectively). Arrows 
along trajectories indicate the forward direction of time. 
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Figure 12: Geometrical relation between 6 zt and for (n, m) = (3, 2.5) and 7 = 10 as in Fig. 11 
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Figure 14: Dependence of the key dynamical variables as a function of the initial condition yo = (yio,0) on the 
yi=axis for (n,m) = (3,2.5) and 7 = 10: exit time i ±fc (yo) on into the non-oscillatory regime beyond the 
intervals Ae^ ±1,=F0 ^ shown in figure 13 (top panel); N turn {yo) (middle panel); and yi c (yo) (bottom panel). In each 
panel, "circle" and "crosses" symbols correspond to points in Ae + ( fc+1 )'~ fe and Ae - ^*" 1 ', respectively. Notice 
the alternate structure in panel c) reflecting the spiralling topology of the boundaries b + and b~ shown in figure 13. In 
order to construct panel c), we have sampled each turn segment Ae ±fc+1,=Ffc by 20 yio points. The two end points are 
chosen to be less than 10~ 8 away from p ±fc + 1 and p Tk . The other 18 points are equally spaced within Ae ±fe+1,Tfe . 
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Figure 15: Scaling laws associated with the fractal properties as a function of initial condition at turn points yo 

m) = (3,2.5) as in Figure 11: a) 7 



into j 0) of b + for yi > on the yi-axis for 
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Figure 16: Comparison of the theoretical predictions (|1 10|), (111), (113), (115) for the exponents a,b,c,d defined 



by (|104j)-(|107|) with estimations obtained from the direct numerical integration of the dynamical equations using a 
fifth-order Runge-Kutta integration scheme with adjustable time step. The lines are the theoretical predictions as a 
function of n for m = 2.5 (solid line), m = 2.75 (dotted line) and m = 3 (dashed line) as function of the exponent n. 
The symbols correspond to the exponents obtained by numerical simulation for 7 = 10 (square), 7 = 100 (diamond) 
and 7 = 1000 (crosses). 
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Figure 17: Non-parametric test of the prediction ( 1 1 8 ) by rewriting ( 1 1 8 ) as ( 1 IS ) so as it qualifies by a linear behavior. 
The time t +k is plotted (triangles) as a function of \y\ k \{t +k ) which are proxies for (A yi (t). The triangles are fitted 
to ( 1 19) to get t* and f3. The exponent b is fixed to its theoretical value given by ( 110|J113|J115| ). The first panel shows 
the whole calculated range. The second panel shows a magnification close to the exit point of the oscillatory regime. 
The different straight lines corresponds to fits of the data with (115) over different intervals, with t* = 0.3857 and 
P = 0.3381 for the interval k = 1 -> 301; t* = 0.6665 and (5 = 0.3381 for the interval k = 61 -► 301; t* = 1.9577 
and (5 = 0.3380 for the interval k = 241 -> 301; t* = 2.4330 and (3 = 0.3381 for the interval k = 271 -> 301; 
i* = 3.5695 and /3 = 0.3379 for the interval k = 286 — > 301. 



56 




0.1 0.2 0.3 0.4 



-50 



||i|iSjji|i \ \ \ \ \ \ 

Mw Ti r/ 7 /' / 



0.1 

^-100 



0.2 



0.3 



0.4 



0.5 



-50 



4 
3 
2 
1 

°\ 
-2 
-3 
-4 



0.07 



0.08 



0.09 



-150 



-100 



\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ N 



:lirrrrrrrrr/r/ 7 7 7 7 7'/ 7 7 7 7 7 7 

-3 
-4 



0.055 



0.06 



0.065 



Figure 18: Terminal critical value y\{t c ) = y\ c as a function of initial value y\§. The top panel shows the oscillations 
of yi c as a function of yio in the first 300 turn segments Ae^ ±( - k+1 ^' Tk ^ (see definition 7.3.9). Recall that y\ c diverges 
at the boundaries corresponding to the intersection of the two curves b^ 1 with the yi-axis. Due to finite numerical and 
graphical resolution, we can only observe a cusp-like behavior associated with points approaching very close to these 
boundaries. The other three panels gives magnifications of the top panel for the first 50 turn segments Ae^ =t( -' c+1 ^ Tfc ^ 
(second panel), from the 50th to the 100th turn segments Ae^ k+1 ^ Tk "> (third panel) and from the 100th to the 150th 
turn segments Ae( ± ( fc+1 ) ,=Ffc ) (fourth panel). In each panel, "circle" and "crosses" symbols correspond to points in 
^ e +(k+l),-k anc j ^ e -k,+(k-i) ^ respectively In order to construct these panels, we have sampled each turn segment 
^ e ±fc+i,q=fc by 20 y w points. The two end points are chosen to be less than 10~ 8 away from p ±fc+1 and p Tfc . The 
other 18 points are equally spaced within Ae ±fc+1,=Ffc . 
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Figure 19: Same as figure i.e., terminal critical value yi(t c ) = y\ c as a function of initial value y\§ = yi(to) 
from the 240th to the 260th turn segments Ae (±(fc+1) ' Tfc) (top panel), from the 260th to the 280th turn segments 
Ae (±(fe+i),TA.) ( m iddle panel) and from the 280th to the 300th turn segments Ae^ k+1 ^ Tk ^ (bottom panel). The two 
vertical lines provide a guide to the eye to verify the almost perfect self-similarity. 
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Figure 20: Same as figure from the Oth to the 20th turn segments Ae^ =t ^ c+1 ^ =Ffc ^ (first top panel), from the 20th to 
the 40th turn segments Ae^** 1 )'^ (second panel), from the 40th to the 60th turn segments Ae^t^ 1 ^^ (third 
panel) and from the 280th to the 300th turn segments AeW^ 1 )' 1 ^ (fourth bottom panel). The two vertical lines 
provide a guide to the eye to show that self-similarity is not qualified. 
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